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A phase-field model that allows for quantitative simulations of low-speed eutectic and peritectic 
solidification under typical experimental conditions is developed. Its cornerstone is a smooth free- 
energy functional, specifically designed so that the stable solutions that connect any two phases 
are completely free of the third phase. For the simplest choice for this functional, the equations of 
motion for each of the two solid-liquid interfaces can be mapped to the standard phase-field model 
of single-phase solidification with its quartic double-well potential. By applying the thin-interface 
asymptotics and by extending the antitrapping current previously developed for this model, all 
spurious thin-interface corrections to the dynamics of the solid-liquid interfaces can be eliminated. 
This means that simulation results become independent of the thickness W of the diffuse interfaces. 
As a consequence, accurate results can be obtained using values of W much larger than the physical 
interface thickness, which yields a tremendous gain in computational power and makes simulations 
for realistic experimental parameters feasible. Convergence of the simulation outcome with decreas- 
ing W is explicitly demonstrated. Furthermore, the results are compared to a boundary-integral 
formulation of the corresponding free-boundary problem. Excellent agreement is found, except in 
the immediate vicinity of bifurcation points, a very sensitive situation where noticeable differences 
arise. These differences reveal that, in contrast to the standard assumptions of the free-boundary 
problem, out of equilibrium the diffuse trijunction region of the phase-field model can (i) slightly 
deviate from Young's law for the contact angles, and (ii) advance in a direction that forms a fi- 
nite angle with the solid-solid interface at each instant. While the deviation (i) extrapolates to 
zero in the limit of vanishing interface thickness, the small angle in (ii) remains roughly constant, 
which indicates that it might be a genuine physical effect, present even for an atomic-scale interface 
thickness. 



I. INTRODUCTION 

The phase-field method has emerged in recent years as a powerful tool to study interface dynamics in areas such 
as solidification , solid-solid phase transformations |3| , and fluid mechanics 0, Q . Their common point is that the 
dynamics of interfaces is coupled to one or several transport fields, which can lead to the spontaneous emergence of 
complex interfacial patterns. The characteristic scale of such patterns is typically macroscopic, whereas the interfaces 
are rough on a scale of a few times the range of the interatomic forces. Because of this intrinsic scale separation, the 
evolution of pattern-forming interfaces has traditionally been formulated in terms of free-boundary problems (FBP), in 
which the interfaces are assimilated to mathematical surfaces without thickness However, the numerical treatment 
of FBPs is highly non-trivial, because one needs to discretize and track the sharp boundaries. 

The phase-field method avoids this need by introducing additional fields to distinguish between phases. These 
"phase fields" take different, constant values in each bulk phase, and interpolate smoothly between these values 
through a diffuse interface of thickness W . The equations of motion for these auxiliary fields are set up so that one of 
their intermediate level sets, which will represent the interface, obeys the desired FBP in the so-called "sharp-interface 
limit", in which W tends to 0. 

To construct these evolution equations, the physics of diffuse interfaces is typically used as a guide, and often some 
qualitative physical meaning can be attached to the phase field. In solidification, for instance, it can be interpreted 
as an order parameter, and its time evolution considered to be a relaxation towards the minimum of a free-energy 
functional, in the spirit of the time-dependent Ginzburg-Landau models for the out-of-equilibrium thermodynamics of 
phase transitions 6]. The sharp-interface limit of the model is then checked a posteriori by matching, order by order, 
terms of asymptotic expansions for the fields in powers of VF in a region of slow (the bulk) and rapid (the interface) 
variation of the fields. At the lowest order in the interface thickness W ^ this procedure is quite straightforward and 
yields a FBP that does not exhibit a dependence on W . 
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However, since the numerical discretization needs to resolve the rapid variation of the fields through the interfaces, 
simulations are only feasible using finite W. Furthermore, in order to reduce the gap between the scales of the 
interface thickness and the pattern, W often needs to be chosen several orders of magnitude thicker than in reality. It 
is therefore mandatory to assess the influence of W on the simulation results. This can be done rigorously by going to 
the next order in W in the above matching procedure. In this "thin-interface limit" , in which W is small but finite, 
an effective FBP is obtained, which now generally contains terms that scale with W. 

While these terms may reflect genuine (but small) physical effects associated with the thickness of the interfaces W, 
they have to be eliminated either by a judicious choice of the model and its parameters 0, or by additional terms in 
the evolution equations tailored to counterbalance these effects Q , or by a combination of both |3| . This is necessary 
to scale up W without altering appreciably the outcome of simulations. If all terms in W can be eliminated, the 
outcome is independent of W, at least below some value of W for which the effective FBP holds (and above which 
terms of higher order in W become important). This complete elimination was achieved for the first time by Karma 
and Rappcl for the solidification of a pure substance with equal thermal diffusivities in both phases (symmetric model) 
0. Their pioneer work made it possible to quantitatively simulate dendritic growth in three dimensions 0- More 
recently, the approach has been extended to the one-sided model, in which the diffusivity vanishes in the solid phase, 
the relevant case for alloy solidification Q . 

Industrial alloys often have more than two components and more than one solid phase. As a first step to deal 
with such alloys, we formulate and test here a phase-field model that carries over the above advances to two-phase 
solidification, a problem that is both practically relevant in metallurgy and of interest for the physics of pattern 
formation: Eutectic and peritectic two-phase composite structures are, together with dendrites, the most common 
solidification microstructures found in industrial alloys |lOj |. At a binary eutectic or peritectic point, two different solid 
phases coexist with the liquid. Under certain conditions, both solids grow simultaneously in a cooperative manner 
(coupled growth). For eutectic alloys, the resulting composite can consist either of alternate lamellae of each solid 
phase or of rods of one phase immerged in a continuous matrix of the other. 

Due to the presence of two different solid phases, coupled growth can exhibit a rich variety of different patterns. 
The best studied setting, used as a relevant example throughout this article, is the directional solidification of thin 
eutectic samples [11. 12. 13tlT^ IT5 | . In this process, samples are pulled from a hot into a cold region with a constant 
velocity, and after a transient the solidification front reaches a steady state characterized by a fairly uniform lamellar 
spacing. This spacing is not intrinsically selected by the system, but depends on the growth history. Upon a (large 
enough) change of pulling velocity, the front undergoes various instabilities and exhibits many classical features of 
nonlinear physics, such as bifurcations to states of lower symmetry, periodic limit cycles, spatiotemporal chaos and 
soliton-like traveling perturbations |0, 0, 0, 0, 0| . 

Numero us p hase-field models for two-phase or multi-phase solidification have been pubhshed [11 El ES IM El 
Ell^, 25, 26, 27]. While such models can qualitatively [25l l2^ l2^ (or, in some cases, even quantitatively for some 
aspects 28J) reproduce experimental observations, more accurate modeling is necessary in order to carry out detailed 
comparisons to experiments on specific substances. Albeit all those previous models have the correct sharp-interface 
behavior, it seems quite difficult to perform a thin-interface analysis of them. As a consequence, the influence of the 
interface thickness on their results remains uncontrolled. 

The reason for this difficulty is intricately linked to the construction of the models: The starting point for any 
asymptotic analysis is the equilibrium front solution that connects two different bulk phases. For a binary alloy, there 
are two or more fields involved (the composition and one or more phase fields) through two or more coupled nonlinear 
partial differential equations, generally without known analytic solution. In the sharp- interface limit W 0, the 
coupling between the equations tends to zero, so that a separate solution of each equation can be found and used 
in the asymptotic analysis; however, in the thin-interface limit {W small but finite), the coupling terms remain. A 
strategy to construct more tractable models is to choose the coupling terms such that they vanish in equilibrium, 
which obviously facilitates finding analytic equilibrium solutions for arbitrary coupling strength. Such a model has 
recently been developed for a dilute binary alloy 0. 

For two-phase solidification, an additional difficulty arises, because one now needs to distinguish between more 
than two phases. The first phase-field models for eutectic growth used the standard solid-liquid phase field, and a 
concentration field 0, 0| or a second phase field [2I to distinguish between the two sohds, a and p. Again, 
the equilibrium front problem takes the form of at least two coupled nonlinear partial differential equations with no 
known analytic solution. Later on, it was proposed to assign a separate phase field pi to each phase i, indicating 
its presence {pi = llor absence {pi = 0). Each phase field can then be understood as a local volume fraction, and 
^■Pi = 1 locally [231 . If, on an interface that connects phases i and j, no other phase is present {pk = 0, k ^ 
the above sum rule permits to eliminate one of the two remaining phase fields from the evolution equations, and 
the interface can be described by a single variable, which, again, makes the existence of an exact analytic solution 
more likely. However, the stable solution for the phase fields across an interface does typically present some amount 
of the other phase(s). Although this recalls the microscopic physical phenomenon of adsorption at an interface. 
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this effect has to be ehminated, because it would scale with the interface width. Recently, the use of the so-called 
double obstacle potential |29| in the volume fraction formulation ,30j has permitted to reduce or even eliminate the 
presence of supplementary phases in the interfaces. In conjunction with the so-called interface-field method this 
methodology has recently been applied to simulate eutectic growth '27*1. However, no thin-interface analysis of such 
models is presently available, and might even be complicated by the singular nature of this free energy. 

Having all this in mind, our strategy is to develop a phase-field model for two-phase solidification with a smooth, 
specifically designed free-energy functional that ensures the absence of third phases in the interfaces and with a 
coupling between the phase and concentration fields that vanishes in equilibrium. Although we use as a starting 
point a free-energy functional, we will not seek to link our model to explicit thermodynamic alloy models. Instead 
of this, the final goal is to simulate the given free-boundary problem (FBP) as efficiently as possible; the free energy 
functional is "engineered" in this perspective. 

We obtain a model that can accommodate arbitrary eutectic or peritectic phase diagrams and that exactly reduces, 
for each of the two solid-liquid interfaces, to the standard phase-field model for the solidification into a single solid 
phase, with the standard quartic double-well potential. Obviously, this is advantageous, since we can build on the 
progress that has recently been made on this model In particular, we extend to two-phase solidification the so- 

called "antitrapping current" , a phenomenological addition to the phase-field equations that was recently introduced 
in order to achieve a quantitative phase-field formulation of alloy solidification 8]; more details will be given in 

Sec lmTTl 

We thoroughly test our model in two-dimensional simulations and compare the outcome to results obtained with 
the boundary integral method • We find that the dynamics of the solid-liquid interfaces are accurately simulated 
and independent of the interface thickness for sufficiently thin interfaces, as predicted by the asymptotic analysis. 
However, the dynamics of trijunction points where the three interfaces meet differs from those traditionally postulated 
in the classic free-boundary formulation. These differences, due to the diffuseness of the trijunctions, are explored in 
detail in the numerical simulations fSec. IIV)l . 

A preliminary account of some of our results has been given in Ref. 32]; here, we give a detailed presentation of 
both the model and the numerical simulations. The rest of the paper is structured as follows: We first recall the 
physics of eutectic solidification and write down the classic FBP (Secnj. Next, we construct our phase-field model 
step by step (Sec. Illlj) : we explore the mapping to single-phase solidification on the solid-liquid interfaces to deduce 
their thin-interface behavior, and then refine the model to make this behavior match as closely as possible the FBP of 
Sec.m In Sec. llVl we present the numerical simulations and the lessons one can learn from them. We then conclude 
with a summary of our main results and a discussion of further prospects fSec. IVj) . 



II. FREE-BOUNDARY PROBLEM 



Consider a binary, eutectic alloy. As a relevant example, we take the transparent organic mixture CBr4-C2Cl6, for 
which accurate and extensive experimental data are available El 111 111 IP ■ Its phase dia gram is shown in Fig. ^ 
where T is the temperature and C the composition, expressed as molar fraction of the solute, C2CI6. The liquidus 
lines of two solid-liquid equilibria meet at the eutectic point (Ce,Te), the lowest melting point of the alloy. At this 
temperature, liquid of composition Ce can coexist with two solids of composition Ca < C'e and Cp > Ce- These 
compositions define the limits of the eutectic plateau, of total width AC — Cp — Ca- 

For slow solidification latent heat rejection can be neglected. The temperature is then constant throughout the 
system for isothermal boundary conditions (isothermal solidification). For directional solidification in an imposed 
thermal gradient, one must further assume equal heat conductivities in all phases for the temperature field to be 
independent of the interface position and fixed by the experimental setup only. For a gradient G directed along the 
z direction with pulling speed Vp in the negative z direction, we then have 

T = TE + Giz~vpt), (2.1) 

where we have chosen the origin of the z axis at the position of the eutectic isotherm a,t t — 0. 
The temperature of a solid-liquid interface is given by the generalized Gibbs- Thomson condition, 

T = T^ + m,{C+-CE)-^^^-- (2.2) 

where C+, k and u„ are the concentration on the liquid side of the interface, the local interface curvature, and its 
normal velocity, respectively, and i — a,(3. Up to the first two terms on its r.h.s.. Equation H2.2|l is the equation of 
the liquidus lines, where nii are the liquidus slopes (toq, < 0, nip > 0), taken at the eutectic point, but assumed to 
be independent of the concentration. The third term on the r.h.s, where are the solid-liquid surface tensions and 
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FIG. 1; Experimental phase diagram of the transparent organic eutectic alloy CBr4-C2Cl6 (modified from J. Mergy, G. Faivre 
and S. Akamatsu). C is the concentration of C2CI6, and we have indicated the eutectic temperature Te, the equilibrium 
concentrations of each phase at the eutectic temperature, Ca, Cp, and Ce, and the width AC of the eutectic plateau. 



Li the latent heats of fusion per unit volume, represents the capillary effect which lowers the melting temperature 
of convex parts of the solid. Finally, the fourth term stands for the effects of interface kinetics; /i; therein is the 
linear kinetic coefficient, defined as the proportionality constant between the local undercooling and the velocity of 
a flat interface. Eliminating the temperature between Eqs. (|2.1|l and 1)2. 2|l . we obtain boundary conditions for the 
concentrations at an interface of specified position, velocity, and curvature [Eq. (|2.5c|l below]. 

Growth is limited by solute transport, much slower than that of heat. For solidification in thin samples, convection 
in the liquid is suppressed, and solute transport occurs by diffusion only. Furthermore, diffusion in the solid is in 
most cases so slow, that the final distribution of impurities is still the trace left by the concentration at the solid side 
of the solidification front C_ , which is assumed to follow the equilibrium partition relation 

C_ = k,C+, (2.3) 

where ki are the partition coefficients. Therefore, we neglect diffusion in the solid phases (one-sided model). Then, 
the velocity of the interface and the diffusion flux on the liquid side are related by the Stefan condition, 

v,,{C+-C-) = -h-{DVC\ + ), (2.4) 

where n is the unit vector normal to the interface and pointing into the liquid. This condition expresses mass 
conservation: The concentration difference between the two phases has to be transported away by diffusion for the 
interface to move. As a corollary, in the one-sided model the solid-solid interfaces cannot move, and hence their shape 
is just the trajectory followed by the trij unctions. 

All of the above already specifies the dynamics of each solid-liquid interface independently. Each of them is thus 
found to obey the following free-boundary problem: 

dtc = DV^c (2.5a) 
Dh-^c\+ = [c, - (1 - fc,)c+]t;„ (2.5b) 

c+ = T ( +d,K + P,vA , (2.5c) 



where we have defined the scaled concentration field 



C-Ce 



AC 

and the scaled equilibrium concentrations of the two solids at the eutectic temperature, 

Ci — Ce 



(2.6) 



(2.7) 
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Equation H2.5a(l holds in the hquid, and Equations l|2.5b() and H2.5c|) on the interface; the minus (plus) sign in Eq. H2.5c(l 
is valid for the a-liquid (/3-liquid) interface, Zi^t is the z position of the interface, 

are the thermal lengths. 



the capillary lengths, and 



(2.9) 



= I ^Ar^ ' (2.10) 

kinetic coefficients. We also define the average capillary and thermal lengths, d — [da + dp)/2 and It = {It + ^t)/-^' 
as well as the diffusion length 

Id = D/vp. (2.11) 

Finally, the a-liquid and /3-liquid interfaces are connected by imposing local equilibrium at trijunction points where 
three interfaces (a-liquid, /9-liquid, and a-f3) meet. The balance of the surface tension at the trijunction points reads 



toLCToL + i/3L0'^L + ta/jCTQ/S = 0, (2-12) 

where tij is a unit vector tangential to the i-j interface and pointing outwards from the trijunction. This condition, 
known as Young's law, yields the equilibrium angles between any two interfaces. An important special case is a steady- 
state composite structure such as the lamellae depicted in Fig. |21 for which the solid-solid interfaces are aligned with 
the z direction. The angles of the i-liquid interfaces with the x direction, called contact angles (we adopt the notation 
of Refs. 13113)) are given by the solution of the equations 

(Jap = (TaL sin 6'a + cr;5L sin 6*^ (2-13) 
CTqLCOS^q = (T/3L cos6'/3. (2-14) 

For equal solid-hquid surface tensions, we have sin^Q, = sinOp = crQ^/(2(TiL)- 

Note that we have not taken into account crystallographic effects, which would appear through anisotropics in the 
kinetic coefficients and the surface tensions; the latter would also lead to additional "Herring torque" terms in Young's 
law. It has been shown by boundary integral simulations |l6l IT^ that all the relevant instabilities and morphologies 
can be reproduced by the above model. Therefore, we have limited ourselves to an isotropic formulation both for the 
FBP and the phase-field model. 

The growth of lamellar composites has been analyzed in the seminal paper by Jackson and Hunt '36] . They found 
that there exists a family of steady-state solutions parameterized by the lamellar spacing A (see Fig.lJJ. The average 
undercooling of the front, AT, depends on the spacing as 

AT=^fA + ^V (2.15) 



where 

ATjH - 



2\/2AC 



, , , , 7p(77)[(1 -??)rfa sin6i„ +77rf/3 sine'^;]/;^ (2.16) 
7^(1 - 77) \ma + mi3 J y 

AjH - ^j2lD[il-v)daSme^+ f]dp sin 9p]/P{v), (2.17) 

with P{r]) — '^'^^iSm^{'Kr]n)/{TTTi)^, and rj the nominal volume fraction of the a phase at the eutectic temperature, 
related to the global sample composition 0^0 by Coo = V^a + (1 — v)'^P 13 HH US- The front undercooling has a 
minimum (AT — ATjh) for the spacing Ajh, which constitutes a reference length for lamellar eutectics. 
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FIG. 2: Sketch of the lameUar geometry with the definitions of the contact angles 9a and 
phase, and A the lameUar spacing. 



rj is the volume fraction of a 



III. PHASE-FIELD FORMULATION 



In this section we construct and justify our phase-field model. We start by introducing our notation and a variational 
formulation fSec. lIIl'^ . Next, we introduce a free-energy functional that provides the minimal features of the model, 
i.e., three distinct phases with interfaces that are free of any adsorbed third phase, and enough freedom to fit a phase 
diagram fSec. IIII'B|) : The equilibrium properties of the model, including the phase diagram, are then verified in Sec. 
nil CI The model is analyzed for a two-phase system, and it is found that it exactly reduces to the usual phase-field 
model of single-phase solidification (Sec. IIII D() . This mapping allows us to deduce the thin- interface behavior of 
two-phase interfaces without performing a detailed asymptotic analysis. Furthermore, this analogy serves us as a 
guideline to extend the model, and, in particular, to allow for two independent kinetic coefficients on the two solid- 
liquid interfaces (Sec. IIII E|) . to introduce a non-variational formulation fSec. IIII F|l . to address the one-sided case in 
which the solid diffusivity is neglected (Sec. IIII G)l . and the case of different surface tensions fSec. IIIIIl)l . 

Although we use as a starting point a free-energy functional, we should stress that we will not seek to link our 
model to an explicit thermodynamic formulation. The final goal is to obtain a model that simulates the free boundary 
problem of the previous section as efficiently as possible; the free-energy functional is "engineered" in this perspective. 



A. Framework and notation 



We use three phase fields pi, i = a, (3, L, and denote p= {pa,Pf3,Ph)- Each field represents the volume fraction of a 
different phase, and hence we would like pi G [0, 1] \fi and 

J2p^^l. (3.1) 

i 

We then introduce a free-energy functional 

F= f fdV, (3.2) 
Jv 

defined as the volume integral of a free-energy density 

/(p, Vp, c, T) = A7grad(Vp) + Hfpip) + XMp, c, T), (3.3) 

whose dependencies have been broken down in the following way: fp depends only on the phase fields and provides 
a "free-energy landscape" in p; in particular, it contains the analog of the double-well potential of single-phase 
solidification, fc couples the phase fields to concentration and temperature and hence defines the phase diagram. 
It does not contain a term in (Vc)^, which turns out to be unnecessary. Both fp and fc are dimensionless; H and 
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X are constants with dimensions of energy per unit volume, /grad sets a free-energy cost for gradients in p, forcing 
interfaces to have a finite width. We wiU use a quadratic form in the gradients of the phase fields, and thus generalize 
the squared-gradient term of single-phase solidification. Therefore, K has dimensions of energy per unit length. 

In a variational formulation of the equations of motion, the non-conserved phase fields are assumed to evolve towards 
the minimum of F, 



..dp, 1 SF 



Vi, (3.4) 



where H has been introduced on the right hand side to remove the dimensions of -F, and t{p) is a relaxation time 
that is the same in all Equations H3.4() (i.e., Vi), but may depend on the local values of the phase fields (see Appendix 
^for a detailed motivation of this choice). The functional derivative has to be evaluated subject to the constraint of 
Eq. 1)3. This can be done by the method of Lagrange multipliers: A term ^(1 — X^i -Pj) added to the free-energy 
functional; the derivatives are taken as if the variables pi were independent, and ^ is determined and eliminated. The 
result for three phases is 

where the variations on the r.h.s. arc now taken as if all pi were independent. We recall that variations are evaluated 
according to 

5F df df 

^S^TJTTf-^, (3.6) 



Spi{x) dpi ^ dipvPiY 

where v counts the spatial coordinates. Equations H3.4|) and H3.5|l automatically ensure Y^idpi/dt = 0, which is 
necessary to maintain consistency with the constraint of Eq. H3.1() . Note that we do not impose at this level pi G [0, 1] Vi, 
in contrast to other formulations 's^. This restriction will instead result from the specific design of our free-energy 
functional. 

The concentration c is a conserved field, and thus obeys the continuity equation 

|+V./=0, (3.7) 

where J is the fiux of scaled concentration c, 

-.5F 

J = -M{p)V — , (3.8) 

with M{p) a mobility. 

As already mentioned in the introduction, the key point of our model is that an interface between any two phases 
i and j should be free of any adsorbed third phase, that is, pk = Vfc 7^ If we construct a free-energy functional 
T such that Pk = 0,1 be the only homogeneous, stationary stable solutions of Eq. H3.4fl for p^, then a i-j interface 
will be free of phase k, since any perturbation around this solution would relax back to zero. More generally, this 
requirement will ensure that bulk phases pk — 1 are stable, i-j interfaces are stabilized only by the presence of bulk 
phases i and j at each side and contain no third phase, triple junctions or lines can only exist where interfaces meet, 
and so on. Note also that, with appropriate initial and boundary conditions, this guarantees that pi e [0, 1] Vi (for 
the continuum equations; in the discretized equations, small overshoots may occasionally appear). The requirements 
on J-^ for stationarity and stability read respectively 



SF 
Spk 
PF_ 



= Vfc, (3.9a) 
> Vfc. (3.9b) 



Pa+Pl3+Pl, = l, Pk=0,l 

Furthermore, in order to ensure stability with respect to concentration fiuctuations, we need to require 

d'^F 

> 0, (3.10) 

which is just a standard condition of thermodynamic stability. 

To construct a suitable free-energy functional, our strategy is to choose the two parts of the free-energy density 
that depend on the gradients only (iiT/grad) and the local values of the fields only {Hfp + Xf^), respectively, such 
that the volume integral of each one separately satisfies the above conditions. 
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B. Minimal model 

For clarity of exposition, we first construct our model taking the parts of the free-energy density that functionally 
depend on the phase fields only, /grad(Vp) and fpif), to be invariant under the exchange of any two phases. We will 
later break this symmetry in different ways as it becomes necessary. 

For the energetic cost of the gradients of the phase fields, /gradj we adopt the simplest scalar expression that is a 
regular and isotropic function of Vp and respects this symmetry, 

/grad = ^5](Vp,;)'. (3.11) 

i 

It is straightforward to check that the volume integral of this function satisfies Eqs. 1)3. 9|l . 
As for the potential part /p, we take 

fp = /tw, (3.12) 

where the triple-well potential /tw is the analog of the double well in the standard phase-field model of solidification 
with a single solid phase. In order to acquire some "geometric intuition" about its construction, it is useful to visualize 
it with the help of the Gibbs simplex. We recall (see also Appendix^ that the Gibbs simplex is an equilateral triangle 
of unit height, where each vertex represents a different phase (either a or /3 solid, or liquid), and, for any point, the 
(signed) distance to the side opposite to a given vertex represents the value of the phase field associated to that vertex. 
Thus, each side corresponds to a purely binary interface, and the center, Pa = Pp — Ph — 1/3, to the triple point. 
Points outside the triangle represent negative values of one or two phase fields. 

The parts of the free-energy density that depend only on the local values of the fields, but not on their gradients, fp 
and /c, can be plotted as a surface over this simplex. Moreover, the conditions of Eqs. H3.9|) on their volume integrals 
can be rewritten using Eq. 1)3. 6|l to replace functional with partial derivatives: Focusing on the part that we are now 
concerned with, fp = /grad, the conditions read dfp/ dpk\p^+pp+pi^=i, pk=o,i = and 9^/p/9p^-|p„+p3+pL=i, p^=o,i > 0. 
It can be easily verified that the left hand side of the first condition is equal to the slope of the fp surface, plotted 
over the Gibbs simplex, in the direction k. Therefore, the conditions of Eqs. H3.9a|l and l)3.9b)l can be understood as 
requiring zero slope (flatness) and convexity, respectively, in the direction normal to a side, both on that side and on 
the vertex opposite to it. In other words, the sides of the triangle must be valleys of the free-energy landscape, so that 
the stable solutions that connect two vertices run exactly along the sides (purely binary interfaces), and each vertex 
should be a local minimum; hence, there will necessarily be a saddle point on each side that will connect the two 
vertices of the side, and each interface will thus "feel" a double-well potential. There should be no other minimum. 

The requirement for minima at pk = 0,1 along the pk direction [Eqs. (|3.9|l ] suggests to construct the triple-well 
potential as a sum of equal double- well potentials /dw(p) for all the phase fields: 

/tw = X! /Dw(pj)- (3-13) 

i 

If fowip) ha-s two equal minima at p = 0, 1, /tw will have equal global minima in each phase i characterized by 
Pj — Sij. Applying the prescription of Eq. H3.5|l for = 0, 1, we find that the slope is equal to {2/2>){dfY)\N / dp){p = 
Pk) - {l/3){dfuw/dp){p = Pi) - {l/3){dfuw/dp){p = Pj). Since pk = 0,1 and {dfDw/dp){p = 0, 1) = 0, the first term 
is zero; furthermore, sd pk = 1 Pi = Pj — and the other two terms vanish too; as for pk ~ 0, then pj = \ — pi, and 
it can be seen that flatness is satisfied whenever /dw(1 ~ p) = fowip), i-e., if /dw is even in 2p — 1. Similarly, it 
can be shown that convexity is also satisfied as long as the double well convexity in the wells overcomes its concavity 
on the maximum in between, 2{d^ fuw / dp^){p = 0, 1) > -~{d^ fBw/dp^){p — 1/2). The standard quartic double well, 
proportional to 

/dw =p'(1-p)', (3.14) 

satisfies the above requirements. This choice is particularly convenient since it will allow us to relate our model to 
the standard phase-field model. The result for /tw is plotted in Fig. O 

So far, all bulk phases correspond to equally deep wells, which means that they all have the same bulk free-energy 
density. However, the relative stability of phases generally depends on the local values of the concentration and the 
temperature. In the phase-field model for single-phase solidification, this is implemented by adding a function to the 
free-energy density that tilts the double well by an amount that is proportional to the local driving force for phase 
transformations. Here, we need to construct a function that tilts the triple well in such a way that the conditions of 
flatness (and, if possible, convexity), are maintained. Note that the condition of flatness in particular, Eq. (|3.9a|l . not 




FIG. 3: Triple well potential, visualized as a surface over the Gibbs simplex. The black lines mark the borders of the Gibbs 
simplex, drawn on the surface; for this function, they also coincide with the trajectory of the equilibrium interface solutions. 



only ensures the absence of third-phase adsorption, but also that the two wells of every pure i-j interface remain at 
Pi = 0, 1 in spite of the tilt, a condition known to be important for single-phase solidification Q. 

We begin by constructing a function gi which raises the well i, gi{pi = 1) = 1, and leaves the other two as well as 
the entire interface between them unaffected, gi{pi = 0) = 0. On the other two interfaces {pj = with j ^ i), we 
require gi to be antisymmetric with respect to the saddle point of /tw at Pi ~ Pj — 1/2, i.e.. 



g^{l~p,,Pi,0) = 1 - g,{p,, 1 -Pi,0). 



(3.15) 



This requirement is necessary to avoid undesirable thin-interface corrections in the matching to the free-boundary 
problem (see Refs. 13, Isll)^ and to adjust surface tensions independently of the phase diagram (see next subsection), 
and is hence important 38]. Replacing 1 —pi by pj and then exchanging the mute indices i and j, we can rewrite the 
requirement in the form 



gj{p^,Pj,0) = 1- g,{pi,pj,0), 



(3.16) 



which can be understood as the analog of pj — 1—pi for pk — 0. This latter form [Eq. (|3.16(l ] will be used below. We 
also impose zero slope on all sides of the Gibbs simplex in accordance with Eq. (|3.9a|l . but not convexity, for reasons 
which will soon become apparent. 

The lowest-order polynomial in the phase fields that satisfies all the above requirements is fifth-order and unique 
at this order, 



5. = ^ {15(1 - p.) [1 - {pk ~Pjf] +P^ (9p2 - 5)} 



(3.17) 



it is plotted in Fig. 0] Remarkably enough, this function reduces to the polynomial pf (10 — 15pi -I- 6p|) on the i-j 
and i-k interfaces when the constraint pi + Pj + Pk — 1 is taken into account, which happens to be proportional to 
the tilting function used in the available quantitative models of single-phase solidification 0, • 

We then couple the phase fields to the concentration c and the temperature T through the functions gi{p)'. We use 
these functions to interpolate between three free-energy parabolic wells in concentration 



for the three phases i = a, /3, L, characterized by pj = gj ~ Siji 



(3.18) 



/c- 2 



-J2MT)g,{p) 



B^(T)gM- 



(3.19) 



The Ai{T) and Bi{T) define the equilibrium phase diagram, as discussed in the next subsection. 
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FIG. 4: Elementary tilting function Qi, visualized over the Gibbs simplex. The i vertex is to the right. 

With fc specified, we can now write down an explicit expression for the concentration current. We define a 
dimensionless chemical potential, 

''=l^'-£ = ^ = '-T.MT)9.ip), (3.20) 

i 

SO that Eq. H3.8|l becomes 

J -D{p)^H, (3.21) 

where D{p) = M{p)X is a phase-dependent difFusivity. Note that the fc in Eq. (|3.19|) or (|3.20|) satisfies the stability 
condition of Eq. (|:ilO(l . 

Equation (|3.9a|l for no third- phase absorption also holds for fc, since each of the gi satisfies this condition by 
construction. However, no general statement can be made about the convexity of fc, because it depends on the form 
of gi and on the values oi Ai, Bi and c. It can be shown that for special choices of the Ai and Bi, the second derivatives 
of fc with respect to the third phase can be made to vanish throughout all borders of the Gibbs simplex. In any case, 
the other components of the free-energy density will ensure global convexity and hence stability for sufficiently small 
values of the ratio X/H. The solutions can always be stabilized by adding "obstacle" terms in the free energy, as 
discussed in Sec. 1111 HI below: however, our experience is that, whenever an instability occurs, the model is anyway 
too far from the thin-interface limit to represent the correct free-boundary problem, so that the convexity of fc is not 
an issue in practice. 

C. Equilibrium solutions 

In order to illustrate how the evolution equations of the model are derived and to clarify the reasons for the particular 
choice of Eq. (|3.19|l for fc, we derive the stationary solution for a planar interface. We check that this solution satisfies 
the conditions of thermodynamic equilibrium, compute its surface tension, and obtain the relationship between the 
coefficients AiiT) and Bi{T) and the cquilibrimn phase diagram. 

1. Evolution equations 

We first write down explicitly the evolution equations for the minimal model, whose building blocks have been 
derived above. For the phase fields, we take the functional derivative of the free energy with respect to each phase 
field Pi according to Eq. 13. 511 . where the component of the free-energy density /grad is given by Eq. I|3.11|l . fp — fxw 
by Eq. (|3.13|) with Eq. 13.14|l . and fc by Eq. (|3.19|l with Eq. H3.17|l . Since c and fi are related, we can eliminate one 
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in order to reduce the number of variables; we choose to eliminate c. We obtain: 



dpi 



9/i 
'dt 



dt ' 



(3.22) 
(3.23) 



where we have defined W = ^/kJH, A = X/H, and 



Po+P/3 +PL = 1 



= [(Pfe-w)'(3p»-2) + (l-p0'(3p» + 2)] , z^j,fc, j^k, 



(3.24) 



9pi 



Pc,+P/3+PL=l 



1 % 

2 ap. 



i+P/3+PL = l 
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(3.25) 



The evolution equation for /i is a diffusion equation with a source term that reflects the redistribution of solute at 
advancing interfaces. 

Next, consider a generic interface. Since we have taken care that our free-energy functional, if not all of its 
components, satisfies Eqs. (|3.9|) for = 0, 1 (fc 7^ hj), the phase field that corresponds to the third phase k vanishes, 
Pk = 0. Therefore, the sum rule of Eq. (|3.1|) can be used to eliminate pj through pj — 1 — pi [and hence gj in terms 
of gi, Eq. (13.161) . gk{pk = 0) = 0], so that we are left with a single phase field pi as desired. Equations H3. 22113. 2^ 
thus yield [33 



t{p} 



dpi 
dt 



W'^'p,-2p,{l-p,){l-2p,) 



- 15Ap2(i - p,)2 (A^. _ Ai) - [B, - Hi)] 



'dt 



^ = V 



dg^jpk = 0) 
dt 



(3.26) 
(3.27) 



where gi {pk ~ 0) is evaluated using pj = \ — pi 



2. Stationary solutions 

We search for stationary planar interface solutions of Eqs. H3.26|l and H3.27|l . that is, solutions that have vanishing 
time derivatives but vary in space along a single coordinate x and connect bulk phases i [pi{x = —00) = 1] and j 
[pi{x — -\-oo) = 0]. From the condition dfi/dt = we obtain that the chemical potential /i, if bounded, has to be a 
constant in space. This is of course one of the conditions of thermodynamic equilibrium, and we hence denote this 
constant by fj^l^^. Its value is fixed by the second requirement, dpi/dt = 0, through a solvability condition for the 
stationary, one-dimensional version of Eq. H3.26|l . 

W^d^^p, = 2k(1 -p,;)(l - 2k) + 15Ap,2(l -k)2 [/i(A, - A,) - {B, - B,)] = (3.28) 

where the second equality defines a "potential" V{pi) = — (1— p^)^ — (A/2)pf (10— 15pi+6pf ) [/i (Aj — Ai) — {B,j — Bi)] 
up to a constant. This notation refers to the well-known "particle-on-a-hill" analogy, in which Eq. (|3.28|) is interpreted 
as an equation of motion with time coordinate x for a point particle of position pi that moves in the potential V . 
Since there is no dissipative term, a solution to this equation that connects pi — and pi — 1 exists if and only if 
V{Q) = V{1), which requires the squared bracket [Aj — Ai) — {Bj — Bi)] to vanish and hence yields 



(3.29) 
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The phase-field profile obtained with the remaining term of Eq. (|3.28() and the given boundary conditions reads 



1 



tanh ( — 1= — - I 
V V2VF J 



the desired solution for an interface centered on . Equations (|3.29l) and 
(|3.2U|) ] and the condition /i = /i^-jj determine the concentration profile, 

c{x) = + A,g,[p{x)] + A^g^[p{x)]. 



(3.30) 

together with the definition of /_* [Eq. 



Taking the limits x ±cxd, we find the concentration of phase i in equilibrium with phase j, 



J-3 — 



A, 



Bj — Bi 
Aj Ai 



(3.31) 



(3.32) 



By choosing appropriate functions Ai{T) and Bi(T), we can reproduce any phase diagram as characterized by its 
coexistence lines (T) . Since only free energy or concentration differences between phases are relevant [note the 
form of the square bracket in Eqs. H3.26|l or (|3.28|) ]. we may choose Aj^ — Bl — without loss of generality. With 
the remaining four functions of T we can fit two liquidus and two solidus lines; the solid-solid equilibrium is then 
constrained, but this is irrelevant for the one-sided case, for which the solid has no dynamics. 



3. Relation to thermodynamics 

To obtain the thermodynamic interpretation of the above solvability condition V{0) = V{1), note that the two terms 
in the middle of Eq. H3.28|l correspond to {l/2){d/dpi)fp{pi,l—pi,0) and {l/2)X{d/dpi)fc{pi,l — Pi,0), respectively 
■ The second is a partial derivative because fc also depends on pi implicitly through c{pi) [Eq. H3.31I1 ]: its 
total derivative is {d/dpi)fc — {d/dpi)fc + {dfc/dc){dc/dpi), where dfc/dc — fi is the constant /i^-'q in equilibrium. 

Therefore, this second term can also be written as {\ /2)X{d/ dpi){fc — jJ-c) . By formal integration of the second equality 
in Eq. H3.28|l we thus obtain (for constant /i) 

V{P^) =-\ {fp{p^. 1 - K, 0) -f A [f,{p,. 1 - P^. 0, c) - ^^ijc] } (3.33) 

up to a constant, to be compared to the (dimensional) grand potential density, 

u^f-X^xc^ X/grad + Hfp + X{f, - ^ic) = A7g,.ad - 2HV{p,), (3.34) 

where the last equality holds for constant /i and we recall that A — X/H. Since /grad vanishes in bulk phases, the 
condition T^(0) = V{1) implies that the grand potential density be equal in the two coexisting phases, uj{Q) ~ '^'(1) = 
'^eq- We hence see that the two conditions for the existence of a stationary solution correspond to equality of chemical 
and grand potentials. Since fp takes the same value in all bulk phases, we note that the equality of grand potentials 
simply requires that of fc — ^J.c. The latter, together with a constant fi = dfc/dc, is what is usually represented 
graphically in the well-known common tangent construction, illustrated in Fig. \^ 



4- Surface tension and choice of fc 

At this point, it is important to realize that the conditions of stationarity, or, equivalently, the common tangent 
construction, only require the grand potential uj and the quantity fc — fJ-c to take equal values at both sides of the 
interface. However, they both generally vary through it, which means that their surface excesses are typically finite. 
The surface excess of the grand potential, for instance, gives the surface tension, 

/ + 00 
{u;[p{x),cix)]-ijg}dx. (3.35) 
-C30 

The particularity of our model can now be understood: At equilibrium, fi = /i*-^ happens to make the square brackets 

and hence the whole term proportional to A in Eq. (|3.28|l vanish for any value of pi, i.e., throughout the whole interface. 
This decouples the phase and chemical potential fields at equilibrium, as emphasized in the introduction. More 
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FIG. 5: The fc part of the free-energy density is an interpolation between three parabohc bulk phase free energies; the 
equilibrium compositions and chemical potentials can be obtained by the common tangent construction illustrated by the 
straight, dashed lines for both solid-liquid equilibria. This construction corresponds to the conditions of equal chemical and 
grand potential in the two coexisting phases. 



precisely, the equilibrium phase-field profile [Ea. (|3.30|l ] balances the derivatives of ii'/giad and Hfp only; moreover, 
since we saw that the term proportional to A corresponds to {l/2)X{d/dpi)fc{pi,l—Pi,0) = {l/2)X{d/dpi){fc — fJ.c), the 
fact that it vanishes means that fc — is constant throughout the interface, i.e., it has no surface excess. Therefore, 
it does not contribute to the surface tension, which again is determined purely by the remaining -fC/grad + H fp. We 
obtain 



a,j = ^^KH = ^ WH. (3.36) 

To better realize the importance of this point, it is useful to consider the roles of the coefficients K, H and 
X in the free energy. X is the dimensional prefactor of the concentration-dependent terms and hence sets the 
magnitude of the thermodynamic driving forces, which is a macroscopic, physically measurable quantity. In contrast, 
K and H can be adjusted in order to achieve a desired surface tension and interface thickness. In our model, 
this is particularly straightforward since Cij does not depend on X: The interface thickness W can be varied for 
computational convenience while keeping the measured value for cr^ just by fitting H [see Eq. (|3.36|l ]. and, still, X 
can be chosen independently to match its measured value. The capillary lengths di are proportional to the ratio of 
the surface tension to the driving forces {di cx an^/X). Since (Ty oc WH, di cx WH/X — W/X, a scaling that will be 
confirmed by the thin-interface analysis given in the next subsection. This means that scaling up simultaneously W 
and A will leave the physics invariant, as desired. 

For a generic model, in contrast, the terms proportional to A do not vanish in the equation of motion for pi nor in 
the expression for the grand potential. Therefore, the phase and chemical potential fields do not decouple; both the 
phase-field profile (and therefore its thickness W) and the surface tension depend on the three constants K, H, and X, 
and no general analytic solution exists. In the expression for the surface tension, there is an extra contribution coming 
from fc, which, on dimensional grounds, should be proportional to WX. Therefore, changing W while keeping the 
surface tension and driving forces fixed is far more complicated. Furthermore, the coupling to fc will introduce some 
dependence of the surface tension and capillary lengths on the composition and temperature. While such dependencies 
may reflect some physical effects, it is difficult to control them properly since they will be blown up as the interface 
thickness W is scaled up to allow for simulations. 

The advantage of our model comes from the fact that the partial derivative of fc with respect to any of the two 
phase fields present on an interface, {dfc/dpi)\p^+pp+p^^i^p^^o, vanishes at equilibrium /or any value ofpi, i.e., thanks 
to the particular fc used. More specifically, we achieve this through an fc whose partial derivative splits in (i) a factor 
which depends on pi only and (ii) one which depends on n only, but none on c. In turn, for the particular type of 
polynomial construction of fc we use, (i) can be traced back to the symmetry of the g^'s [Eq. (|3.16|l ]. and (ii) to the 
fact that the second derivatives of the bulk-phase free energies with respect to concentration are constant and all 
equal. The symmetry of the gi^s can be understood physically as ensuring that fc depends only on relative differences 
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of free energy and concentration between any pair of phases i and j, for any phase fraction of the two phases pi and 
1 — Pi. This is reflected in the expression in square brackets in Eqs. H3.26|l or H3.28|l . Apart from this symmetry, the 
gi are arbitrary, phenomenological phase-field interpolation functions. 

In contrast, the free energies of the different bulk phases are, in principle, functions that are fixed for a given alloy 
system and cannot be freely adjusted. Our free energy is therefore an approximation, justified by the fact that it 
yields the desired phase diagram. However, apart from the phase diagram, the free energy also determines the latent 
heats, as we shall outline below, and these enter the capillary lengths [Eq. H2.9|l ]: moreover, due to the Gibbs-Thomson 
effect or to kinetics, the concentrations at both sides of a solid-liquid interface can deviate from the prediction of the 
phase diagram at a given temperature, which affects the amount of impurity actually rejected. As we shall see, our 
approximation for constraints the capillary lengths and slightly modifies the impurity redistribution; however, the 
first is compatible with the experimentally measured values, and the second is completely negligible for low-speed 
solidification. To get rid of this approximation, it should be possible to use an approach recently introduced for dilute 
alloys in which internal energy and entropy are interpolated by two different functions However, the generalization 
of this method to two-phase solidification, that is, the introduction of two different sets of gi functions, is outside the 
scope of the present paper. 



D. Mapping of each solid— liquid interface to single-phase solidification 

Here, we address the behavior of our model for small but finite values of the interface thickness W. It would be very 
interesting to investigate this behavior including the triple junctions; however, the problem there is quite involved, 
since two independent phase fields and the chemical potential vary rapidly and are all coupled. Therefore, we limit 
ourselves to an analysis of the interfaces, taking advantage of a mapping to single-phase solidification. The behavior 
of the trijunction points will be investigated numerically and discussed in detail in Sec. IIVI 

We are actually interested in solid-liquid interfaces, since the solid-solid one has no dynamics. We hence set 
j = L in Eqs. (|3.2()|) and 13.27|l . To proceed, we rewrite the expression ^(^l — Ai) — {B\^ — Bi) thus appearing 
in Eq. (|3.26() as (^l — — Meq)j where = (Bl — -Bi)/(^L — ^i) is the dimensionless chemical potential for 

coexistence of the solid i and the liquid as given by Eq. (|3.29|) ]. It hence becomes apparent that the driving force, 
— lh\pj{l — piY{Ai^ — Ai)(^yi — /Zpq), is proportional to a deviation from equilibrium. Let us first consider isothermal 
solidification, so that Ai(T), Ai^i^) and /i*q(r) are just constants. Then, the evolution equations can we rewritten in 
terms of a new variable 

(/^-/4q)/(^L-A,), (3.37) 

which measures departure from equilibrium on the i- L interface, and a new phase field (^i = 2pi ~ 1, which takes the 
values +1 and —1 in the i solid and liquid phases, respectively. We obtain: 

= W^'V^^. + 0,(1 - 0?) - ^(Al - A.if\{l - <\>\fu (3.38) 



du -± 



i^(p)V« (3.39) 



where j denotes the other (absent) solid phase. 

With a constant relaxation time, r(p) = tq, and a constant diffusivity, D(j)) = D, the above equations constitute 
precisely the phase-field model for the solidification of a pure substance treated by Karma and Rappel, if the com- 
bination (15/8) (^L — Ai)"^}, is identified with the coupling constant A used in Ref. 0; more precisely, we recover 
the variational version of their model. Therefore, we can use their results on its thin- interface behavior, i.e., the 
behavior for values of W much smaller than a typical lengthscale of the microstructural pattern. This is described by 
an effective free boundary problem for the field u. Undoing the change of variables of Eq. H3.37|l and applying Eq. 
(13.20(1 . we can compare it with the desired free-boundary problem of Eqs. 1(2. 5|l : 

We find that the Gibbs-Thomson boundary condition, Eq. ((2. Set , is satisfied, with capillary lengths and kinetic 
coefficients related to the phase-field parameters by 

W 

\Ai^~Ai\X 



Pi = ai 



a2\Ai^ - A, 



(3.41) 



where ai = -\/2/3 and 02 = 0.7464 are numerical constants related, but not equal to ^3 those in Ref. 0. We 
note that, in general, the expressions for the capillary lengths and the kinetic coefhcients are not the same for the 
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two solid-liquid interfaces, but depend on the A^'s, and hence on concentration differences [recall Eq. (|3.19|) ]. As 
anticipated, the parameter A controls the ratio of the interface thickness W to the smallest physical lengthscales, the 
capillary lengths di. The former can be varied while keeping the latter fixed just by changing A accordingly. For 
quantitative output, the simulation results should hence become independent of A for A small enough. This will be 
checked for our model in Sec. IIVI Note that, as discussed in detail in Refs. 0,13, convergence may be achieved for A 
much larger than 1, since it is only assumed that W is much smaller than a typical length scale of the pattern, which, 
in turn, is usually much larger than di. 

The diffusion equation, Eq. H2.5a|l . is also satisfied. As for the mass conservation condition, we find 



Dn-{Vc\+-Vc\.)^{A~Ai^) 



(3.42) 



where Vc|+ and Vc|„ denote the concentration gradients on the liquid (+) and the solid (— ) side of the interface, 
respectively. There are two differences with the Stefan condition of the desired FBP, Eq. (|2.5b|) . Obviously, in the 
one-sided case the diffusion flux on the solid side vanishes, and hence the second term on the left hand side is absent. 
We address this case in Sec. IIII Gl But there is a second difference: The right hand side of Eq. (|3.42|) corresponds 
to the right hand side of Eq. I|2.5b|l only for constant concentration gaps, ki — 1. To see this, substitute Eq. (|2.2|l . 
absorbing its temperature dependence in the phase-diagram concentration gaps cf^ — c}^, into Eq. Ij2.5b|l . This yields 
Dh-Vc = Vn [c*^ - ± (1 - ki){din + /3iU„)] . The term c*^ - c% corresponds to our Ai - Al [see Eq. (|S22l], but 
the other terms are lacking. This is due to the use of equal fc,i/dc^ for all three bulk-phase free-energy functions 
fc^i [Eq. H3.18|l ]. as announced in Sec. IIII CI When the chemical potential is shifted from its temperature-dependent 
equilibrium value for a flat interface, the corresponding shifts in concentration are determined by d^/dc = fddc^ 
and are hence the same at both sides of each solid-liquid interface. This explains why our model displays the right 
deviations of c from its phase-diagram value at the liquid side [correct Gibbs-Thomson condition, Eq. (|2.5c|l ]. but 
does not reflect them in the concentration gaps that appear in the impurity rejection. However, the magnitude of the 
missing terms is small, since the departure from the equilibrium phase diagram is very small in slow solidification: 
The PiVn term is usually neglected and will be set to zero in our simulations; as for diK, a typical curvature k is 
given by the inverse of the lengthscale of the microstructural pattern. In turn, this lengthscale goes as y/dfo, where 
Id = D/Vp is the diffusion length. Therefore, diK ~ ^fdjTo. For typical experimental values of slow solidification 
djlo ^ 10^"' — lO^'"^, so this correction can be safely neglected. 

The other constraint coming from the equal fc.il dc^ is reflected in Eq. I|3.40|l : The ratio of the capillary lengths 
is fixed to 



da _ \Ap-Ai,\ 
dp 



\Ac 



At 



(3.43) 



This relation can be understood from thermodynamic considerations: The latent heats of the two solid-liquid phase 
transformation can be evaluated from the bulk-phase free energies fi of [Eq. (|3.18(l ] by Li = r(sL — Sj), where 
Si = —{dfi/dT)\c is the specific entropy of phase i (and similarly for the liquid). Making use of the conditions of 



equal chemical and grand potentials, we find Li — {diJ'^/dT)\ 



solid-Hquid equihbria (dfi^^/dc) = d'^fc/dc'^\cq are the same, we find La/Lp 
Using the definition of the capillary lengths, Eq. (|2.9|l . this yields 



(a/i^y9c)|cf - c}^\/mi. Since for both 

{cT^-cl^)/icf~C)\{mp/m^). 



da 




„/3L 




_ CqL 


Ap-Ai^ 


dp 










Aa-Ai^ 



(3.44) 



where the second equality makes use of Eq. I|3.32|l . Since, in our minimal model, ctql — c^l, Eq. (|3.43|) follows. In 
general, the ratio depends on the temperature through the {T). Near the eutectic point where solidification occurs 
for small undercoolings. 



da_ 

dp 



crph 



(3.45) 

to within experimental 



This latter relation is satisfied by the reference alloy of our present study, CBr4~C2Cl6, 
accuracy. 

Note that these lacking terms in the solute rejection and the restriction on the ratio of the capillary lengths do 
not vary with the interface thickness, but are zeroth-order corrections: They are not introduced by the diffuseness of 
the interfaces, but by our approximation for the free energy, which can be replaced by more sophisticated choices if 
needed, as already discussed at the end of Sec. IIII (T4l 
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Let us conclude this section with a comment on directional solidification, where the assumption of constant Ai{T) 
and Ai^{T) is no longer valid. If we still want to apply the mapping to single-phase solidification, we must neglect 
the temperature variation of Ai{T) and Ai^{T) within the region where the mapping holds. The scale of variation of 
the temperature is set by the thermal lengths defined in Sec. ^ whereas the mapping region needs to extend over 
at least a few times W in order for the thin-interface behavior of the model to be well defined (and the results of 
Ref. to apply). Because these results already assume that the lengthscale of the pattern is much larger than W, 
and taking into account that the thermal lengths are typically far larger than this pattern lengthscale, the condition 
^ is automatically met. Therefore, even in presence of a temperature gradient, the results of Eqs. H3.4(J|1 and 
(I3.41() apply, with the value of T (and hence Ai{T) and j4l(T)) taken in the center of the interface. 



E. Independent kinetic coefficients 

The relation between the phase-field parameters tq, A and W and the interface kinetic coefficients of the FBP is 
given by Eq. H3.41|l . where \Ai^ — Ai\ — — \cf^\ are set by the phase diagram, _D is a material parameter, and ai 
and a2 are numerical constants. Given a desired accuracy as fixed by W or, equivalently, A [linked by Eq. 1)3.40(1 ] . 
the only free parameter left is the phase- field relaxation time tq. However, we need at least two free parameters in 
the model, since there are two kinetic coefficients that are a pr^or^ independent. In the formulation of Ref. (H) 
with a double-obstacle potential, this is implemented by specifying different relaxation times for each binary {i- 
j) phase transformation. In our formulation with a smooth free energy functional, this is not possible, as shown 
in Appendix ^ Instead, we achieve independent kinetic coefficients by keeping the same relaxation time r for all 
possible transformations, hut making r dependent on the phase fields. 

In order to maintain the mapping to the model of Ref. Q on the individual interfaces, r has to be a constant along 
those interfaces. A function taking different, constant values on all three interfaces necessarily has discontinuity lines 
inside the Gibbs simplex. Here, we are mainly interested in controlling the dynamics of the solid-liquid interfaces, 
since, in the final one-sided model, the solid-solid interface does not move. Therefore, we use a function that takes 
constant values (tq and r^) on the two solid-liquid interfaces only, and interpolates smoothly between them inside 
the Gibbs simplex: 




^13 



t{p) = { ' ^ 2 (p„+p,) if PL ^ 1 ^3 

if PL = 1, 

where f = (tq, -|- Tp)/2. Note that the only discontinuity is on the vertex pi^ = 1. It did not cause any problems in 
practice, but, if needed, it could be smoothed out in a small neighborhood of the vertex without inducing appreciable 
errors in the calculations. 

As a result of this procedure, the two kinetic coefficients. 



Pi = Oi 



W 

a2\Ai^-A; 



]Ai^-A,\\W 'D 
can now be adjusted independently by choosing the two constants and 



(3.47) 



F. Non- variational model 

It has been shown for phase- field models of single-phase solidification '7,'3?| that it is often advantageous to introduce 
an additional degree of freedom by switching to a non- variational formulation, which exploits the two different roles of 
the tilting functions gi: In the evolution equations for the phase fields, Eq. (|3.22l) . they favor one bulk phase over the 
other; in the one for the chemical potential /i, Eq. H3.23|l . they constitute a source or sink of /i, which corresponds to 
impurity rejection or adsorption, respectively. The gi need to satisfy certain common requirements for both of their 
roles (although for different reasons), namely to interpolate from to 1 as pi goes from to 1 and to be antisymmetric 
with respect to the point pi — pj = 1/2 on i-j interfaces [Eq. (|3.15|l ]. However, the requirement of flatness, Eq. 
(|3.9a() . constrains only the derivatives with respect to the phase fields, which do not enter the evolution equation of 
the chemical potential fj,. 

Therefore, we can switch to a different chemical-potential-like variable. 



^ = c - ^ Aj/i,(p), 

i 



(3.48) 
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which amounts to replace the gi in Eq. (|3.20|) by new functions hi that have the same hmits and symmetries as the gi, 
but do not necessarily satisfy Eq. H3.9a|l . This implies that the dtgi in Eq. (|3.23ll are replaced by dthi. Furthermore, 
a different expression for /c has to be used to derive the phase-field equations, namely 

/c = E 5^ (p1 (^) ' (^)] ' (3.49) 

i 

which is the form given in our preliminary account js^ . The mapping to single-phase solidification can be repeated: 
For the choice hi = Pi, the result corresponds exactly to the isothermal variational model of Ref. 0; as a result, 
Eqs. (|3.4Q(I and H3.47|l hold again, with the same value of ai as before, but now with 02 = 1.175. 

The choice hi(p) = pi is quite advantageous for simulations, because the source term in Eq. (|3.23|) is less peaked 
inside the interface than for the fifth-order polynomial gi{p), which makes it possible to use a coarser discretization, 
resulting in a considerable gain in computational speed, as discussed in detail in Ref. (7(| . 

G. One-sided model 

We now address the one-sided case, for which the solute diffusivity in the solids is neglected compared to that of 
the liquid D. This is much more realistic for alloy solidification than the symmetric model considered up to now. We 
hence set 

D{p)=Dq{pi^), (3.50) 

with a function q{ph) that interpolates between 1 in the liquid and in the solid. 

This seemingly small change has important consequences. We already mentioned in Sec.^that a strictly vanishing 
solid diffusivity prevents solid-solid interfaces to move. For the solid-liquid interfaces, the introduction of a phase- 
dependent diffusivity breaks the symmetry of the impurity diffusion between solid and liquid. Furthermore, the 
variation of the diffusivity D{p) through the diffuse solid-liquid interfaces needs to be considered in the asymptotic 
analysis, which becomes quite involved ^.37J; we will only resume the most important points here. 

It turns out that, in general, the effective free-boundary problem (FBP) obtained from the thin- interface analysis 
displays several terms that are not present in the original one and that scale with the interface thickness W. As 
explained in detail in Ref. , in analogy with the thermodynamics of diffuse interfaces these terms can be linked to 
"surface excesses" , which are the integral through the interface of the excess of a quantity, where excess means the 
difference of its actual value at a point inside the diffuse interface region and its reference bulk value at the nearest 
side of the interface (the bulk values in the solid and the liquid might differ) . 

Let us consider the evolution equation for the chemical potential on the solid-liquid interfaces, Eq. 13.27|l . where the 
gi were replaced by hi as explained in the previous subsection. Three relevant surface excesses have been identified. 
The first is that of the "source function" hi{p): since it determines the equilibrium solute profile, an excess in this 
quantity confers a net impurity content to the interface (solute adsorption at the interface) . The quantity of adsorbed 
solute increases with the area of the interface, which modifies the mass conservation condition, Eq. (|2.5b|) . at moving 
curved interfaces (interface stretching). For the one-sided model, we have two more surface excesses: that of diffusivity 
D{p}, which leads to solute diffusion along the interface (surface diffusion) and thus also modifies the mass conservation 
condition; and that of the ratio hi(p)/ D{p), which generates the well-known solute-trapping effect and shows up in 
the modified FBP in the form of a jump of the chemical potential through the interface 

For the symmetric model in which D{p) is actually a constant, the choice of hi{p) odd with respect to the center 
of the interface ensures both the absence of interface stretching and of solute trapping. However, for the one-sided 
case, one has to eliminate surface diffusion too. The easiest way to do this is to also choose D{p) odd. The simplest 
possible interpolation, q(j)h) = PL, already satisfies this condition, and will be adopted in the following. The problem, 
however, is then that the surface excess of hi{p)/D{p) does not vanish, so that solute trapping is present. 

With only two adjustable interpolation functions, one can generally not eliminate all three undesired surface excesses; 
more freedom is needed in the formulation of the model. One way to obtain this is to add a term counterbalancing 
one of the three effects, while choosing hi{p) and D(p) to eliminate the other two. Here, we generalize to two-phase 
solidification one possible such term, originally introduced for single-phase solidification in Ref. [8|: We add to the 
solute current a new phcnomcnological contribution, an "antitrapping current" which drives the otherwise trapped 
solute into the liquid phase. 

By analogy with Ref. Q, its form on each solid-liquid interface is easy to determine: It is directed parallel to the 
outward normal of the solid-liquid interface, and its magnitude is proportional to the interface thickness W and the 
intensity of the solute release which caused the trapping, (^l — Ai)dhi/dt. How to interpolate this current between 
the two solid-liquid interfaces is less clear, since no asymptotic analysis or mapping to single-phase solidification is 
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available near a trijunction. We numerically tested various possibilities for interpolations, and chose the one that 
yielded the best convergence properties. It can be heuristically motivated as follows: The outward normal to each 
phase is given by hi — — Vpi/|Vpi|. Let us consider fia- On the a-liquid interface, it is antiparallel to riL, since 
Pa = ^ — Pl; in contrast, on the a-f3 interface it is antiparallel to fip, since Pa = 1 — p/j. Upon growth, the a front 
advances along the direction of ha with a speed proportional to dpa/dt. The antitrapping current should always be 
directed towards the liquid, even if the solid formed is a "mixture" of both solid phases as happens in a trijunction; 
therefore, we choose to direct it along hj^ throughout the whole system. However, only the component of the growth 
speed directed toward the liquid should contribute to the magnitude of the antitrapping current; therefore, we multiply 
the source strength for each phase by the scalar product —til • hi, which is equal to 1 on the solid-liquid interfaces, 
but smaller than 1 inside the trijunction. 
The total concentration current thus reads 

J = -Dpiy^i + 2aWhi^ (^L - Ai){-hi^ ■ fii)^, (3.51) 

where 2a is a prefactor to be adjusted, and we have taken hi = Pi- Indeed, the particular form of the anti-trapping 
current given here is to be used only in conjunction with the simplest choices both for the diffusivity, as given in Eq. 
(|3.50() with q{ph) = Ph, said the source function, hi = pi] see Ref. ^\ for further details. 
With this new concentration current [Eq. H3.51|l ] , we rederive Eq. (|3.39|) : 



du ^-L /I — 0i -4 \ 1 d(j)i / d(hi ^<hi 



dt \ 2 2 dt \ dt 



(3.52) 



We recover the quantitative phase-field model for one-sided solidification with a constant concentration gap given in 
Ref. 13 . Again, this has the advantage that we do not need to analyze the small- 14^ behavior. From Ref. [g we learn 
that the results of Eqs. (|3.4UI) and (|3.47l) for di and /3j and the discussion thereafter still apply, and that a — 1/(2a/2) 
is the value for which the term exactly counterbalances solute trapping. Interestingly, the values of ai and 02 stay also 
the same as given before, for reasons explained in Ref. @. The fact that the only quantitative model for one-sided 
solidification available so far uses hi = pi and that this choice enables one to recover the same numerical values for 
ai and a2 as in the case of constant diffusivity is obviously another reason, on top of its lower numerical cost, to use 
hi = Pi instead of hi — gi . 



H. Unequal surface tensions 

As we have seen in Sec. IIII C 41 thanks to the particular form of our coupling fc between phase and concentration 
fields, the surface tension, or equilibrium surface excess of the grand potential / — pc, reduces to the surface excess 
of the remaining two terms: /giad and fp. Since so far both /giad and fp = /tw are symmetric with respect to the 
exchange of any two phases, the surface tensions of all interfaces are equal. In order to treat unequal surface tensions, 
this symmetry needs to be broken by modifying either of these two terms. In multi-phase-field models, different 
gradient energies are used for each interface (3Q |. In our approach, we want to maintain at least the condition 
of flatness, Eq. (|3.9a|) . and it turns out to be easier to modify the potential part (see Appendix IbI for a detailed 
discussion). 

We add a new term /saddle to this potential part, 

fp = /tw + /saddle, (3.53) 



/saddle = 2^ /saddle, i- (3.54) 
i 

/saddle is duc to chaugc the height of the saddle points on each binary interface by a tunable, different amount, 
through the elementary functions /saddle, which, in turn, should shift the saddle of /tw separating phases j and k 
(i ^ J, k). Therefore, /saddle, i should vanish on interfaces other than the j-k one, and respect the valley character of 
all interfaces. 

A function positive everywhere that vanishes on all sides automatically satisfies the latter two requirements, but 
obviously does not shift any saddle. However, it will be useful later on. The simplest such function is 



/obs =pIpIpI, 



(3.55) 
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FIG. 6: Function used to lift the saddle between two phases and hence alter the surface tension. It is flat along the two 
unaltered interfaces. The position and strength of the hill in the middle are controlled by the parameter b in Eq. (13.5611 (here, 
h= 12). 

which corresponds to an elevation (obstacle) on the triple point and outside the Gibbs simplex. By dropping the 
p1 factor we get a function which actually raises the saddle on the j-k interface and still vanishes and has valleys 
on the others, but which is not flat on the j-k interface in the direction perpendicular to it. We hence make 
the ansatz /saddle, i — P jPfc /saddle, and impose flatness at pi ~ 0. The resulting condition for /saddle, i reads 

3(9/saddie/i9pi)|p„+p^+PL=i = 2/saddie/(PiPfe) at = 0, SO that the simplest choice turns out to be /saddie,i = 
'^PjPk + Spi. This actually corresponds to a function that is small everywhere but in the neighborhood of the j-k 
saddle, where it is maximal. Therefore, its j-k side is flat but concave. To correct this concavity, we add to it the 
obstacle function /obs given above, 



through the term in b. This composite function now has a valley that also runs along the j-k side as long as & > 9/2. 
The larger b, the closer the maximum of this function is to the triple point and the further to the j-k saddle. In Fig. 
|S1 we show this function for 6 = 12. 

On any purely binary interface pi = 0, the whole function /saddle reduces to 2aip^{l — p)^, where p is any of the 
other phase fields. If we now put together the triple-well potential and the saddle functions, we find a free-energy 
density fp = /tw + /saddle = 2p'^(l - p)^[l + aip{l - p)] on such interfaces. 

It is important to note that this modification does not affect the coupling between phase fields and concentration, 
and hence all the equilibrium compositions and the chemical potential remain the same as in the minimal model. In 
contrast, the surface tension is modified. Making use of the equipartition relation (i.e., the fact that /grad = /tw for 
an equilibrium interface solution) , we find the total surface excess per unit area of /grad + /tw to be 



Of course, for = the result reduces to Eq. (|3.86() . For a,; ^ 0, it is straightforward to evaluate this integral 
numerically and to tune the different surface tensions by using different a^. However, the equilibrium profile of a 
j-k interface whose saddle has been shifted by a finite Ui is not any more the usual hyperbolic tangent solution of 
Eq. , but a kink-shaped profile that has to be calculated numerically. Since the equilibrium profile is the starting 

point for the entire asymptotic analysis, all calculations, and in particular the determination of the numerical constants 
oi and a2 that appear in Eqs. (|3.4U|) and 13.47|l . have to be repeated for this new profile, taking also into account the 
extra terms that are generated by /saddle in the equations of motion. Similarly, the form of the antitrapping current 
needs to be adapted, using the methods of Ref. 0. While this procedure is, in principle, straightforward, it is outside 
the scope of the present paper. 

At this point, it is important to realize that the two solid-liquid surface tensions in many eutectic alloy systems 
are quite similar, and only the solid-solid surface tension is markedly different. This is, for instance, the case in 
CBr4-C2Cl6. Therefore, using equal a-liquid and /^-liquid surface tensions is usually a good approximation; only the 
solid-solid surface tension needs to be modified, and hence we have Oq, = = 0, but ^ 0. Then, we recover the 



saddle, z 



aiP]Pk{'2'PoPk + iPi + bp}), 



(3.56) 




(3.57) 
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quantitative phase-field model on both solid-liquid interfaces; the solid-solid interface is not a problem, because it 
does not move (except close to the trijunction) in the one-sided model. 

We should mention here that the computational effort to simulate the situation with unequal surface tensions 
(ofc > for any fc) is generally larger than for equal surface tensions (at = Vfc). This is due to the fact that the 
free-energy landscape of a i-j interface whose saddle has been risen (a^ > 0, k ^ is steeper than that of an 
unaltered interface, so that a smaller grid spacing Ax is necessary to discretize it in a convergent manner. 



IV. NUMERICAL TESTS 



A. Implementation 



Our goal here is to show how the model derived in the preceding section is used in practice and to assess its validity 
and precision. We start by writing down the complete evolution equations for the one-sided, non-variational model 
with antitrapping current. They are essentially the same as Eqs. 13.2211 and (|3.23l) of the minimal model, except that 
Eq. I|3.22|l now includes a contribution from /saddle, D(p) = Dpi,, the gi have been replaced by hi = pi in Eq. (|3.23|l 
but not in Eq. (|3.22(l , and the total concentration current is now given by Eq. (|3.51(l , which includes the antitrapping 

term. We make the equations dimensionless by scaling lengths by T4^ (a; = x/W, z — z/W , V — WV) and times by 
the f = (tq + ti3)/2 used in the equation that sets t{p), Eq. (|3.46(l (< = t/f). The result reads 



dpi 
di 



V^P^ 



-2p,(l - 2k) + J2pjil-pj){l - 2pj) 



+ aiPjPk [{Pj + Pk)fa,i - 2/6,,;] + ^ ajPiPk [{pi - 2pk)faJ + /hj] 



(4.1) 



dpi 



(4.2) 



where fa.k = PzPj +Pk{^ + bpk/3), fbM = PiPj(l/2 + bpk/3), and we recall that the (95j/(^Pj)lpc+Pfj+PL=i are given by 
Eqs. (|3.24(l and (|3.25|l . that ni — Vpi/\Vpi\, and that we use a — l/(2\/2) to exactly counterbalance solute trapping. 

Only dimensionless parameters remain in the equations, namely, the already dimensionless coefficients Ai, Bi, a^, 
6, and A, and the newly defined D = Df/W"^ and f{p) = T{p)/f [with its Hmiting values = Ta/f and Tp = r^/f], 
where t{p) is given by Eq. H3.46|l . The coefficients Ai{T) and Bi{T) are just constants for isothermal solidification, 
but depend on space and time (through the temperature) for directional solidification. Their explicit expression, given 
when choosing the phase diagram in next subsection, will contain the rescaled pulling speed Vp — Vpf/W and thermal 
lengths l{ = If^/W. 

We numerically integrate Eqs. 1)4. l|l for p^ and Pj3 only (the equation for pl is redundant) as well as Eq. (|4.2|l : 
PL is eliminated everywhere through pL — ^ — Pa — Pfi- We use a simple Euler, forward-time, centered-space finite- 
difference scheme with a grid spacing Ax and a time step At slightly below the stability limits of both the discretized 
diffusion and the phase-field equations, the lowest of which reads Af = (l/4)(Ax/Vr)^ min{l/i3, Tq, f^}. We use 
standard second-order accurate finite differences, except for the Laplacians, which are discretized by a nine-point 
formula involving nearest and next nearest neighbors to reduce lattice anisotropy. We use Ax = 0.8 unless otherwise 
stated, which is the largest value for which most results are numerically converged. Exceptionally, some parameter 
sets require Ai — 0.4. 

We simulate two-dimensional directional (or, in one case, isothermal) solidification in a rectangular simulation box 
of total size Ux x Uz, where Ux and Uz are the number of grid points in the direction perpendicular and parallel 
to the thermal gradient, respectively. We consider perfectly periodic lamellar arrays. A minimal representation of 
this geometry consists of a simulation box with no-flux boundary conditions in both the x and z directions and two 
adjacent half lamellae, one of each phase, with their centers on the box boundaries and the trijunction in the middle. 
Typically, we start with completely flat interfaces, and the phase fields are initialized as step functions, located at 



21 




FIG. 7: Illustrative surface plot of the two phase fields pa and pp, for two half lamellae as described in the text. 40 x 40 grid 
points shown; the system is actually larger in the z direction. 

some initial guess for the interface position. These step functions then quickly relax to the smooth solutions for the 
phase fields, while the interfaces begin to curve and drift to adjust their average undercooling. Alternatively, the 
outcome of a previous simulation may also be taken as initial condition. An example for the resulting configuration 
of the phase fields is shown in Fig. [7| It can be clearly seen that the fields smoothly approach their bulk equilibrium 
values outside the interfaces, and that, on each i-liquid interface (at x =const. sufficiently far from the solid-solid 
interface), the other solid phase, j ^ i, is absent (pj = 0) through the entire interface (as z varies). 

In order to reduce the computational effort, we take advantage of two circumstances. First, as can be seen in 
Fig. Q most grid points correspond to bulk phases, so that the evolution equations can be replaced there by simpler 
versions, which leads to an enormous saving in the total number of operations per time step. Furthermore, most of 
the remaining grid points correspond to purely binary interfaces, so that all the terms in the third, vanishing phase 
field can be dropped. Only in the small area near a trijunction do the full equations need to be integrated. According 
to the conditions Eqs. H3.9I) . wherever a phase field takes a locally constant value of zero (purely binary interface 
or bulk phase) or unity (bulk phase), it will remain constant, so that it does not need to be updated, and its whole 
evolution equation can be dropped. In practice, we first compute the Laplacian of each phase field everywhere, and 
then proceed with further calculations for a phase field pi only if the modulus of its Laplacian exceeds a certain 
(small) threshold. Similarly, the source and antitrapping terms in the evolution equation for the chemical potential in 
Eq. (|4.2I) are evaluated for the locally varying phase field(s) only. In particular, this means that they can be dropped 
for bulk phases. Therefore, in the bulk only the simple diffusion equation needs to be solved; in the one-sided model, 
no field needs to be updated in the solid, be it bulk or a — (3 interface. 

Second, since the diffusion and thermal lengths are much larger than the lamellar spacing, the number Uz of grid 
points required in the direction of the thermal gradient is typically much larger than ■ A simple analytical calculation 
(see for example Ref. |36j) shows that the lateral gradients of the diffusion field decay with the distance to the growth 
front over a length of the order of the interfacial pattern (~ Ux), whereas in the direction of the pulling the solutal 
boundary layer decays only on the scale of the diffusion length. Therefore, the spatial resolution can be decreased with 
the distance to the growth front beyond a scale comparable to Ux- This can be achieved by multi-scale algorithms. 
We have adapted the random walker algorithm of Ref. to eutectic solidification; however, a disadvantage of this 
method is its numerical noise. For our relatively simple geometry, it is straightforward to implement a finite-difference 
scheme that uses coarser and coarser grids away from the interface. We used a hierarchy of in total three grids of 
increasing length in the z direction as sketched in Fig.|Sl The solutions on the different grids are connected by simple 
linear interpolation. At the end of the coarsest grid, lateral concentration gradients are negligibly small; therefore, 
we solve the diffusion equation in one dimension beyond this point and connect the solution to the coarsest grid. We 
carefully checked that the dependence of the solution on the position of the interpolation boundaries is negligible. 
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FIG. 8: Configuration of the multi-grid scheme. In front of the solids that grow upward, a fixed number n/ of rows is treated 
on the fine grid with spacing Aa;, which corresponds to a distance d = UfAx. Two coarser grids of spacing 2Ax and AAx 
have again Uf rows each, which corresponds to distances 2d and 4d, respectively. Beyond the coarsest grids, the solution is 
one- dimensional . 



Finally, we also follow the solidifying front by advancing the simulation box whenever the isotherms have advanced 
by once the coarsest grid spacing. Since, far enough in the solid, the diffusivity vanishes and no evolution takes place, 
the coldest part of the solid can be removed without altering the simulation; new points are added on the liquid side. 
The boundaries between grids are also adjusted, and new initial values on the border of the finer grids are obtained 
by interpolation from the coarser grids. 



B. Choice of parameters 



The parameters that characterize a given physical situation can be grouped into different classes. First, there 
are the characteristics of the phase diagram and some materials parameters that depend only on the alloy system. 
Second, there are the control parameters accessible to the experimentalist, namely the sample composition and either 
the undercooling for isothermal solidification, or the pulling speed Vp and the temperature gradient G for directional 
solidification. Finally, there is the lamellar spacing, which cannot be directly controlled in experiments, but can be 
fixed in the simulations by choosing appropriate initial and boundary conditions. 

For coupled eutectic growth at low solidification speeds, we can make several simplifying assumptions and approx- 
imations in order to keep the number of relevant parameters to a minimum. Since the average temperature of the 
solidification front is very close to Te, we can (i) use a phase diagram linearized around the eutectic point, as specified 
in Eq. (|2.2|) of the FBP. This can be implemented by setting the parameters Ai and Bi to 

A, = 0, A=c.T ^^^-'^^l'^-^ = T {h ~ (4.3) 

Ai(T~Tv) , Z — Vr,t 

i3L = 0, B, = T ^ ^tA^-^, 4.4 
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TABLE I: Materials parameters of the transparent alloy CBr4-C2Cl6, as measured in Ref. 33], and the values used for the 
present study (see text). The di listed under Ref. |^] were obtained using the basic definitions of Eq. 12.9|l together with the 
thermophysical data from that reference. The value of Ce used corresponds to rc = 2.5. 



Quantity 


Symbol 


Ref. [22J 


Value used 


Liquidus slope of a 


ma 


-(81 ± 5) K/mol 


-82 K/mol 


Liquidus slope of /3 




(165 ±5) K/mol 


164 K/mol 


Composition of a at Te 




(8.8 ± 0.4) mol% 


8.8 mol% 


Composition of [3 at Te 


Cr 


(18.5 ± 0.9) mol% 


18.5 mol% 


Eutectic composition 


Ce 


(11.6 ± 0.6) mol% 


11.571 mol% 


Capillary length of a 


da 


(9.5 ± 2.0) nm 


9.29 nm 


Capillary length of j3 


d/B 


(3.5 ± 1.0) nm 


3.71 nm 


Contact angle of a 


6a 


(70± 10)° 


30°-66° 


Contact angle of f3 




(67± 10)° 


30° -66° 


Partition coefficient of a 


ka 


0.75 


1.0 


Partition coefficient of (3 




1.5 


1.0 


Impurity diffusivity 


D 


(0.5 ± 0.1) X 10"^ mVs 


0.5 X 10-^ mVs 



where the second equalities are obtained using Eq. (|2.1|) for the temperature and the definitions of the thermal lengths, 
Eq. H2.8|l . Furthermore, we can (ii) neglect the relative variation of the concentration gaps with respect to their values 
at Te- This is equivalent to set ki = 1 (parallel liquidus and solidus lines) in Eq. Ij2.5b|) [4^ . The advantage here 
is that two parameters (the ki) are eliminated and that the Ai become independent of temperature. Finally, at low 
solidification speeds we can also (iii) neglect the kinetic undercooling of the growth front, i.e., we adopt 13a — — 0, 
which is also the most reasonable choice with the information available, since these coefficients are unknown. All of 
these approximations were also made in the classic Jackson- Hunt analysis (3^. 

Furthermore, we use equal surface tensions for both solid-liquid interfaces (ctql = o'/3l), which implies equal contact 
angles da = Op = 9. This is a reasonable assumption with the available data for CBr4-C2Cl6 9a = 70° ± 10° and 
9p — 67° ± 10° [33 ■ In our test calculations below, we use values of the computational parameters for which angles 
ranging from 9 = 30° (all surface tensions equal) to 66° (case of CBr4-C2Cl6) are expected. 

With approximations (i) and (ii) above, the phase diagram can be characterized by the two ratios 

Tni = mp/nia, (4.5) 
rc = \cp/ca\. (4.6) 

Note that we have cp — Ca — 1 due to the normalization, and hence it follows that cp — rd (rc-l-l) and Ca = — l/(rc + l). 
The measured properties for CBr4-C2Cl6 are listed in Table According to the constraint of Eq. H3.45II . for equal 
surface tensions we should have da/dp = rc- From TableHl we find da/ dp = 2.7 ± 0.5, while Tc = 2.5 ± 0.5. Thus, 
the constraint is respected to within the accuracy of the measured data. For our simulations, we have to choose a 
set of values that exactly satisfies the constraint and is compatible with the experimental error bars. We use the 
concentrations to fix Vc — da/ dp = cp/ca = 2.5, and take as a reference length the average 

d={da + dp)/2 (4.7) 

of the measured capillary lengths, d = 6.5nm. The capillary lengths obtained using this value of d together with 
Tc = 2.5 and the constraint are listed on the last column of TableQ] Similarly, we choose = 2 and fix a corresponding 
pair of values for the liquidus slopes within the error bars. Finally, we adopt the mean measured value of the solute 
diffusivity, D = 0.5 x 10~^ m^/s. The measured and used values of the parameters are summarized in Table 

The control parameters which can be varied in an experiment enter our simulations as follows. The global sample 
composition is fixed by imposing the appropriate value for the chemical potential ^(z — )■ 00) = Cqo (recall that = 0) 
as a boundary condition on the liquid side of the simulation box. For isothermal solidification, the undercooling enters 
through Eqs. (|4.3I) and (|4.4|1 in dimensionless form, — (Te ~T)/{miAC). For directional solidification, the pulling 
speed Vp and the thermal gradient G enter through the same equations via the diffusion and thermal lengths Id = D/vp 
and = niiAC/G. Note that ifj^/li^ — r^, so that — [2/(1 + r„i)]lT and l!^ = [2rm/(l -I- rm)]lT, where we recall 
that It = {It + It)/'^ is the average thermal length. Finally, the lamellar spacing A (not to be confused with the 
coupling constant A) is fixed by the lateral size of the simulation box. 

To summarize, the physical conditions to simulate are completely specified by the two asymmetry parameters Vc 
and rm and four lengths, namely, d. It, Id and A. From these four length scales, three dimensionless parameters can 
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TABLE 11: Simulation parameters used in the runs for the symmetric model alloy. Physical parameters are lo/d = 51200, 
It /Id = 4, A = Ajh, and the sample is at eutectic composition. Space and time units are W and r, respectively. 

X/W b \ Vp It/W At/f 

32 10.633 36.197 0.0079729 5334.5 0.012038 
64 5.3164 18.098 0.0019932 10669 0.024077 
96 3.5442 12.066 0.00088588 16003 0.036115 
128 2.6582 9.0491 0.00049831 21338 0.048153 



be constructed. We choose lu/d, A/ Ajh, and It /Id- The second choice is motivated by the fact that the Jackson- Hunt 
minimum undercooling spacing is a reference length for eutectic pattern formation. Dividing both sides of Eq. (|2.17() 
by d yields Ajh/^ oc y^lo/d, where the proportionality constant depends only on the sample composition (through 
the volume fraction rj) and on (through the ratios di/d). Therefore, specifying A/Ajh and lo/d fixes implicitly X/d. 

All physical conditions fixed, we now choose the only truly free computational parameter, the interface thickness 
W. The relevant scale of the pattern is of course the lamellar spacing A, and therefore the resolution of the phase-field 
simulations is given by the ratio X/W. Let us first outline the procedure to determine the parameters when the 
spacing is given in physical units (meters). Then, X/W directly fixes W in meters. Next, the coupling constant A is 
determined from Eqs. (|3.40() for the capillary lengths in terms of phase- field parameters (recall that Ai^ = 0), 

A = ^i^if^ + ^V (4.8) 



d 2 \\Aa\ \Ap^^ 

Finally, the relaxation times (in seconds) are obtained (for arbitrary kinetics) using Eq. H3.47|l . 

r. = A|A.r(^ + a.^), (4.9) 

which yields f = (tq + Tp)/2 and the fi — r^/f. The diffusion coefficient D, pulling speed Vp, and thermal lengths 
needed in Eqs. H4.1|l and (|4.2I) are then obtained by scaling the corresponding dimensional quantities by W and f. 

For isothermal solidification, the dimensionless undercooling is directly obtained from the actual temperature T 

and the alloy properties Te, rm and AC. 

An alternative way to obtain the simulation parameters is to start directly from the dimensionless ratios. Indeed, 

specifying X/W for fixed X/d fixes directly W/d and hence A from Eq. (|4.8|) . For Pa = I3p = 0, Eq. (|4.9|) yields 

Tp/fa = (I A^|/|^a|)^. For the case ki — 1 considered here, Ai = Ci, and hence |A;3|/|^q,| = Tc, and we obtain directly 

fa =2/(1-1- r^) and t/s = 2r^/(l -t- r^). With these values fixed, Eq. H4.9|l with the Ci given after Eqs. (|4.5|l yields the 

scaled diffusivity, 

D = (l/2)a2A(l + r2)/(l + r,)\ (4.10) 

Since W/d a.nd lo/d aie both known, the scaled pulling speed can then be inferred from Id/W = D/vp] finally, the 
scaled thermal lengths are obtained via the ratio It /I d for directional solidification, or the dimensionless undercooling 
Ai is directly plugged into the AiS and BiS for isothermal solidification. 

To illustrate the above procedure, we show in Tables UTI and ITTll the computational parameters for different values of 
the interface thickness for two series of test simulations that will be discussed in detail below. Both sets of simulations 
are carried out at the eutectic composition and for a spacing A — Ajh • The first is for a model alloy with a symmetric 
phase diagram, = Tc = 1. For such an alloy at the eutectic composition, all equations are completely symmetric 
with respect to the interchange of the two solid phases. In this case, we use It /Id — 4 and lo/d = 51200. The 
second set of simulations is for the phase diagram of CBr4-C2Cl6, rc — 2.5 and = 2, and we use lo/d — 41796 
and It/Id — 4, which corresponds to a pulling speed of Vp ~ 1.8 /im/s and a temperature gradient of G ~ llOii'/cm, 
both fairly typical experimental values. 

Finally, the Oi are dimensionless parameters that control the ratio of surface tensions through the relative heights 
of the free-energy barriers between bulk phases. They hence provide a handle on the contact angles. The choice of 
the ai is independent of all the other parameters. For equal solid-liquid surface tensions, we recall that = = 0; 
ol will then be varied in the next subsection to tune the ratio of solid-liquid to solid-solid surface tensions, and thus 
change the solid-liquid contact angle 9. A change from ol = to cl = 12 is necessary to span contact angles from 0° 
to 66° as desired. On the other hand, b needs to be tuned to ensure the convexity of the free-energy landscape; we 
test values ranging from 6 = 3to& = 12. 
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TABLE III: Simulation parameters used in the runs for CBr4-C2Cl6. Physical parameters are lo/d — 41796, It/Id ~ 4, 
A = A.jH, and the sample is at eutectic composition. Space and time units are W and f, respectively. 

\/W D X Vp l^/W 4/W At/f 

32 14.852 42.713 0.013141 3013.7 6027.4 0.0086186 
64 7.4258 21.357 0.0032854 6027.4 12055 0.017237 
96 4.9505 14.238 0.0014602 9041.1 18082 0.025856 
128 3.7129 10.678 0.00082134 12055 24110 0.034474 

C. Isothermal solidification and contact angles 

We begin our simulations by testing whether our model reproduces the correct contact angles at the trijunction 
point. For that purpose, we need to extract first the interface shapes, then the trijunction position, and finally the 
angles between interfaces from the simulations. We proceed as follows. 

The difference pj — pi is computed everywhere, for every pair of phases i ^ j ■ Wherever one of these combinations 
changes sign, we interpolate the position of a i-j interface point from the two adjacent values of Pj — Pi] we also 
interpolate the value of the third phase pk at that point. If pk < 1/3, the point lies on a "true" i-j interface; 
otherwise, it is located inside the third bulk phase on the "prolongation" of the interface beyond the trijunction point. 
In most cases, the latter points are not plotted, and the i-j interface terminates at the trijunction; however, as we 
will see below, it is sometimes useful to plot the prolongation as well, because it can yield information about the 
internal structure of the trijunction. The average undercooling of each interface is obtained from the z position of all 
of its points. In order to avoid sharp cutoffs and discretization errors at trijunctions, the contribution of each point 
on the Pi — Pj isocontour to the undercooling is weighted by 1 —pk, rather than simply counting points on the "true" 
interface and discarding points on the prolongation. 

When all three interfaces enter and leave a particular elementary grid square exactly once, their entrance and 
exit positions computed by the method described above define one straight segment for each interface. If all three 
segments intersect with each of the other two, a triple point is considered to be detected at the average position of 
the three intersections, provided that this lies within the grid square. While this procedure yields an excellent subgrid 
resolution for the position of the trijunction points, the determination of the angles from the slopes of the segments is 
subject to large grid effects (that is, the values obtained are found to depend on the position of the trijunction with 
respect to the grid points). Therefore, we use a more precise procedure: Around a detected trijunction point, a few 
interface points on both side of the trijunction are first obtained with high precision using nonlinear interpolations. 
Then, two such points on each side of the trijunction are used together with a nonlinear interpolation to obtain the 
angles at the trijunction position. With this procedure, the uncertainty on the angles due to grid effects is reduced 
to about 0.1°. 

Before engaging in directional solidification, we test both the model and the above procedures by simulations of 
isothermal solidification, in which the temperature is set to a constant value below the eutectic temperature; in the 
free-boundary problem, the terms (zjnt — Vpt)/l^ are then replaced by constants — where = (Te — T)/{miAC) 
is the dimensionless undercooling of phase i {i — a, (3). Since the interface temperature is fixed, the quantity to 
be selected is now the interface velocity. This velocity is negative for low undercooling (the lamellae melt), positive 
for high undercooling, and exactly zero for a critical undercooling that depends on the lamellar spacing through the 
curvature of the solid-liquid interfaces. For the symmetric model alloy at the eutectic composition, an exact solution 
for this critical steady state is known: Since the composition in the liquid and the temperature are uniform and the 
velocity is zero, the interface curvature is constant [see, for instance, Eq. (|2.2|l ]. so that the two solid-liquid interfaces 
form circular arcs that intersect at the trijunction point. For a given contact angle and lamellar spacing, the critical 
undercooling is 

Ac = 4sin6lY- (4.11) 
A 

We conducted series of simulations with different undercoolings and extracted the velocity of the solid-liquid 
interface once it had reached a steady state; the critical undercooling was obtained by seeking the zero crossing of the 
velocity. The grid spacing was fixed to Ax = 0.4 to better resolve the neighborhood of the critical undercooling (small 
velocities), and the parameters of the first line in TablelTll fexceot Vp and It) were used. We tested five different values 
of cl for which equilibrium angles from 9 = 30° to 6* = 66° are expected from Young's law, with a constant h = 12. 
The predicted angle was then plugged into Eq. H4.11|l to obtain the theoretical value of the critical undercooling. 
The latter is then compared to its measured value, as summarized in Table HVl We find excellent agreement between 
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TABLE IV: Critical undercooling for stationary lamellar states for various values of al (and hence of the contact angle 
6 = 12 in all cases, and the other parameters are given in the first line of Table HH 





9 expected Ac 


(simulation) Ac(6') (theory) 


error 





30° 


0.001624 


0.001628 


0.25% 
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39° 


0.002047 


0.002055 


0.39% 
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48° 


0.002390 


0.002408 


0.75% 
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56° 


0.002686 


0.002713 


1.00% 


12 


66° 


0.002946 


0.002987 


1.37% 
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FIG. 9: Lamellar states, simulated with a constant temperature that yields almost zero growth speed, and comparison to the 
analytical critical stationary solution that consists of circular arcs. Two different values of b are used together with ol ~ 12 
(predicted angle 9 ~ 66°). Since, at constant temperature, the origin of the z axis is arbitrary, the solutions have been shifted 
in order for lamella tips to coincide. The curves shown are the isocontours pa ~ pl and pp = pl, to be compared with the 
theoretical circle arcs with the predicted contact angle 9 (solid curve). The isocontour pa = P/3, exactly vertical and located at 
x/W = 8, has been omitted for clarity. 



simulations and theory; the error in the undercooHng increases with aL (and hence with the contact angle), but 
remains of the order of 1% even for ol = 12, which corresponds to 6* = 66°. In particular, this implies that Young's 
law is satisfied. 

However, if we measure the contact angles directly using the numerical procedure outlined above, we find a consistent 
result only for cl = 0, that is, the measured contact angles are always very close to 30°; for ol ^ 0, we find angles that 
differ widely from the expected values. Furthermore, whereas the detected critical undercooling is almost independent 
of the parameter b, the measured angles for cl = 12 vary by 20° when b is varied between 3 and 12. 

The reason for this behavior becomes apparent in Fig. El We show the analytical solution, which consists of two 
circular arcs with contact angles 6 — 66°, as well as the interfaces extracted from the simulations for two different 
values of b. Whereas the agreement between the two simulations and the theory is excellent far from the trijunction, 
the three shapes start to become different at a distance of about W from the trijunction. The lines pi — pj exhibit 
curvatures which differ from the curvature of the circular arcs of the analytical solution, and which depend on the 
value of b] the position of the trijunction point also depends on b. This behavior is due to the changes in the potential 
energy landscape induced by the variations of b, which affects the internal structure of the trijunction region. Note, 
however, that the radius of curvature of the solid-liquid interfaces far from the trijunction does not change. This 
implies that the critical undercooling remains the same, and hence that Young's condition and the contact angles 
extrapolated from the interfaces to the trijunction are correctly implemented in all cases. 

The lesson from all this is that the local procedure outlined above does in general not yield the "macroscopic" 
angles, but values that are influenced by the internal structure of the trijunction region. A different possibility to 
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obtain the "macroscopic" angles would be to construct the intersection of the two solid-liquid interfaces, extrapolated 
from their shape outside the trijunction. However, this procedure is prone to large errors out of equilibrium, since the 
curvature varies along the interfaces. As a consequence, contact angles are very hard to measure in out-of-equilibrium 
situations when 6 ^ 30°. When 9 = 30°, the "local" procedure works properly. The particularity of this value is, of 
course, the symmetry between all three phases, which entails that the prolongation of the interface inside the third 
phase runs exactly in the middle between the two other interfaces and is not "deflected" from the direction in which 
it enters the trijunction. In contrast, for ol 7^ 0, the local shape and curvature of this isoline depend on the internal 
structure of the trijunction, which explains the errors made in the "local" measurement of the angles. 

Since we are interested in precise information on the angles, we restrict all the following simulations to the case 
of equal surface tensions, where reasonably accurate measurements of the angles can be carried out. An additional 
advantage of this choice is that it avoids the reduced grid spacing needed to resolve the steeper interfaces associated 
to a higher surface tension. For ol = 12 (which corresponds to the actual value of 9 in CBr4-C2Cl6), Ax needs to be 
divided by a factor of 2 (Ax = 0.4) with respect to ol = 0, which implies using a 4 times larger system size (in terms 
of grid points) running for 4 times more time steps (recall that Ai c>c Ax^), which hence takes 16 times more CPU 
time. Having tested the equilibrium angles, in the remaining subsections we focus on out-of-equilibrium simulations 
of directional solidification. 



D. Steady state lamellae 

The convergence of our model to the thin-interface predictions is tested by performing series of simulations for 
fixed physical parameters and decreasing interface width. For each simulation, the interface profiles and interface 
undercooling are monitored to check when the steady state is reached. For comparison, the same series of runs is also 
performed without the antitrapping current, and for an earlier version of our model presented in Ref. that uses 
different interpolation functions gi. Since those functions are not antisymmetric with respect to the point pi =_Pj = 1/2 
on i-j interfaces [i.e., they do not satisfy Eq. H3.15|l ]. this model exhibits several thin-interface corrections [Mls^. 

The parameters used for the simulations of the model alloy with symmetric phase diagram are fisted in Table |TT| 
Simulation times on a 2.4 GHz Intel Xeon processor range from half an hour for the lowest resolution {X/W = 32) 
to three days for the highest {\/W — 128). The results for the interface profile are shown in Fig. ^1 together with 
results of boundary integral calculations performed with the code of Ref. JJi]. Whereas for the present model all 
curves for X/W > 64 superimpose perfectly, for the other models large errors appear. It can be seen that they are 
smaller for the model that only lacks the antitrapping current; however, a decisive progress is made only when all 
thin-interface corrections are eliminated. Even so, a small difference with the boundary integral prediction remains 
in the W/X limit. However, a close examination reveals that the two solutions are simply shifted with respect to 
each other, with is due either to a residual interface kinetics in the phase- field model or to the approximations of the 
boundary integral method. The relative error in the average undercooling is about 0.3 %. Since the advantage of the 
complete model over the other formulations is obvious, it has been used exclusively for all the remaining simulations. 

We also examine the contact angles at the trijunction point for the same series of runs with the complete model. 
We find small deviations from the expected equilibrium value of 30°. In Fig. ^2 we plot the difference, 30° — 9, as 
a function of W/X. It can be seen that it extrapolates to zero in the W/X — > limit. This indicates that there is a 
finite-interface-thickness correction to the non- equilibrium contact angles that vanishes in the sharp-interface limit. 
Since the difference is smaller than 2° for X/W > 64, this effect is invisible in the profiles and does not appreciably 
influence the results for the undercooling. In particular, it cannot be responsible for the 0.1% undercooling mismatch 
described above, since the latter does not extrapolate to zero in the W/X — > limit. 

We now consider the dependence of the average undercooling of the front AT on the lamellar spacing A, which is 
known as the Jackson-Hunt curve, Eq. H2.15|l . We perform a series of runs for the same materials and computational 
parameters as above, but with varying box size (and hence lamellar spacing). For the resolution, we take Xjn/W — 32 
and 64, since, on the basis of Fig. Iior al. we expect results to be converged for Xjn/W > 64. The curves thus obtained 
are plotted in Fig. ^1 for the flrst (Ajh/W^ — 32, squares) and second {Xjh/W = 64, circles) set of parameters of 
Table ITTl The undercoofing AT, scaled by toAC, can be directly obtained through AT / {mAC) = (zint)/^T (note that 
for the symmetric model alloy, the two liquidus slopes and hence also the two thermal lengths are equal, ma = rn(3 = m 
and — l!^ = It), where (zmt) is the z position of the solid- liquid interface, averaged over the lateral coordinate x. 
The line shows the best fit of the higher resolution data to the Jackson-Hunt law, Eq. (|2.15|l . with ATjh and Ajh as free 
parameters. The fit yields Ajh/^d — 0.02468 and ATjh/(»tiAC) — 0.003023, whereas the theoretical values calculated 
from Eqs. (|2.17|) and (|2.16|l are 0.02403 and 0.003251, respectively; this corresponds to relative differences of 2.7 % 
and 7.0 %, respectively. This good agreement is especially noteworthy because in phase-field models which exhibit 
thin-interface corrections the difference between simulated and calculated minimum undercooling spacings is usually 
much larger |1^ ,2^ 22^ ■ Note that the Jackson-Hunt prediction for the minimal undercooling and its corresponding 
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FIG. 10: Convergence test for lamellar shapes in the symmetric model alloy at the eutectic composition. The solid-liquid 
interfaces are plotted for (a) the complete model, (b) without the antitrapping current, and (c) a model with several thin- 
interface corrections taken from Ref. |2fl |. The parameters for all runs are given in Table HTl Thin solid lines: \/W = 32, dotted 
lines: A/VK = 64, dashed lines: \/W = 96, dash-dotted lines: \/W — 128; thick solid line: result of the boundary integral code 
of Ref. 0. 
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FIG. 11: Difference between the contact angle and the equilibrium value 30° as a function of W/X for the same series of runs 
as Fig. [TUT a). 
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FIG. 12: Dimensionless average undercooling ol the Iront versus dimensionless spacing {Id = D/vp is the diffusion length) for 
the first two parameter sets of Table ITTl The curve is a best fit to the Jackson-Hunt law (see text for details). 




FIG. 13: Convergence test for lamellar shapes in CBr4-C2Cl6. The parameters of all runs are given in Table Uni Thin sohd 
fines: X/W = 32, dotted lines: X/W = 64, dashed lines: X/W = 96, dash-dotted fines: X/W = 128; thick solid line: resufi of 
the boundary integral code of Ref . [l^ . 



spacing is based on an approximate description of the front. We know from Fig. llUf a) that the difference in average 
undercoohng between the phase-field model and the boundary-integral method is only 0.3 % near the minimum 
undercooling. Therefore, the 7 % difference with the Jackson-Hunt value must be attributed to the approximations 
in the theory, and not to the phase-field model. 

Next, we repeat a similar procedure for the phase diagram of CBr4-C2Cl6. The parameters are listed in Table. ITTll 
and the resulting interface shapes shown in Fig. 1131 The behavior is qualitatively similar, but the convergence is 
slower. It can be seen that successive interface shapes fall closer and closer together as the resolution is increased. For 
the larger a lamella, the shape is converged for X/W > 64; for the smaller /3 lamella, tiny differences remain visible 
even between 96 and 128. This is not surprising: convergence is harder to achieve when smaller physical details (such 
as the small (3 lamella here) must be resolved. 

Again, the phase-field result is fairly close to the boundary integral calculation. However, this time a clearly visible 
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FIG. 14: Trijunction rotation angle (f) (see text for definition) versus W/\ for the same series of runs as Fig. 1131 

difference remains: Whereas tlie pliase-field result is above the boundary integral for the a lamella, it is below for the 
P lamella. Therefore, the difference is not a simple shift between the two solutions. An explanation for this fact comes 
from the detailed examination of the angles at the trijunction point. Contrary to the completely symmetric situation 
where a and P are equivalent, the solid-solid interface does not remain straight up to the trijunction. Rather, when 
the trijunction is approached this interface starts to tilt and exhibits a small angle with the vertical axis. The 
angles between interfaces remain consistent with Young's law, but the whole trijunction zone is slightly rotated. Since 
we are in a steady state, the trijunction necessarily moves vertically (in the lab frame) as solidification proceeds. 
This implies that the solid-solid interface is not the mere trijunction trajectory, but must move after the trijunction 
passage, i.e., at its solid side, which is possible in a phase-field model, since the diffusivity decays to zero in the 
solid over a distance of order W. In the sharp-interface formulation, in contrast, where the diffusivity sharply falls 
to zero in the solid, the solid-solid interface cannot move and has to be strictly vertical up to the trijunction. This 
difference between phase-field model and sharp-interface description induces the observed shape difference between 
the converged phase-field solution and the boundary- integral result. We plot in Fig. the rotation angle (/) as a 
function of W/X. Contrary to the deviation of the contact angle from its equilibrium value found in Fig.^] this global 
trijunction rotation angle does not extrapolate to zero for W/X 0, but rather seems to approach a finite value. This 
behavior will be discussed in more detail below. 



E. Oscillatory limit cycles 

As it is well known, lamellar steady-state growth becomes unstable above a critical spacing that depends on the 
growth conditions and the alloy composition. At the eutectic composition, the first instability to occur is a period- 
preserving oscillatory instability, that is, all lamellae of the same phase start to oscillate in phase. This instability can 
be captured by our reduced simulation geometry, because the center lines of each lamella remain symmetry planes of 
the pattern, whereas the time translation symmetry of steady-state growth is broken. The oscillations amplify in the 
course of time and finally saturate and reach a limit cycle with a well-defined amplitude. The bifurcation diagram 
of final oscillation amplitude versus lamellar spacing was found to be slightly supercritical in Ref. . Simulating 
oscillatory limit cycles close to the instability threshold therefore constitutes a particularly sensitive test of our model, 
because even slight deviations in the lamellar spacing or differences in behavior will induce large changes in the 
oscillation amplitude. 

We used the symmetric model alloy with the same physical parameters as in the preceding subsection, except the 
lamellar spacing, which was fixed to to A = 2.2 Ajh. The simulations were started with a slightly asymmetric initial 
condition in which one of the lamellae was slightly ahead of the other. The limit cycle rapidly emerged; an example 
for the simulated structure once it has reached its final oscillation amplitude is shown in Fig. 115b for X/W — 64. We 
define the amplitude A of the limit cycle as the maximum lateral deviation of the trijunction from its steady-state 
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FIG. 15: Oscillatory limit cycle in the symmetric model alloy for A = 2.2 Ajh. (a) Snapshots of the two solid-solid interfaces 
(thick lines) and successive solid-liquid interfaces (thin lines), (b): Oscillation amplitude A/\ vs inverse interface thickness, 
\/W . (c): Blowup of the trijunction region. The arrow in (c) indicates the direction of the trijunction motion. 

value, normalized by the lamellar spacing, A — Axmax/A. We plot A versus the resolution \/W in Fig. 115b. It can 
be seen that the convergence is slow: the results still depend on the resolution even for \/W — 192. 

One possible reason for this slow convergence can be recognized in Fig. 115b . where we plot both a snapshot of a 
trijunction region (solid curves) and the final solid-solid interface left behind (dashed curve). Like for the asymmetric 
steady-state solutions, the trijunction seems to be slightly rotated, that is, the direction of the instantaneous solid- 
solid interface is not parallel to the direction of motion of the trijunction point (indicated in Fig. 115b by an arrow and 
roughly parallel to the final solid-solid interface). 

Figure seems to imply that our phase-field model commits large errors when calculating limit cycles. However, 
as already mentioned, we have chosen here a particularly sensitive test case. In Fig. 1161 we show the whole bifurcation 
diagram, that is, the amplitude of saturated limit cycles versus the lamellar spacing, calculated with two different 
resolutions. As shown in the insets, the two curves can be superimposed quite nicely by a simple shift along the A 
axis. Therefore, the relatively large difference in oscillation amplitudes comes from a small shift in the instability 
threshold, while the overall shape of the bifurcation diagram is quite robust. The relative error committed on the 
instability threshold is about 3 %. Hence, our model remains a useful tool for the exploration of nonlinear behavior 
in eutectic solidification, even if the trijunction behavior prevents completely converged calculations. 

F. Discussion 

The most interesting outcome of the above numerical tests is the information about the behavior of trijunction 
points. As already mentioned, to our knowledge no detailed rigorous treatment of the equilibrium or dynamics 
of a diffuse trijunction is currently available, and it seems quite challenging indeed. Therefore, only a qualitative 
explanation of our numerical findings can be given, namely by comparing the phase-field and sharp-interface pictures. 

In the sharp-interface model. Young's law is assumed to hold even outside of equilibrium; furthermore, the solid- 
solid interface cannot move, since the solute diffusivity vanishes in the solid immediately adjacent to the trijunction. 
In the phase-field simulations, both of these assumptions are relaxed. As for Young's law, there is a deviation from 
the equilibrium contact angles. In our simulations of the symmetric model alloy, this deviation decreases with the 
interface thickness and extrapolates to zero in the (numerical) sharp-interface limit. A likely origin of this effect is 
the relaxational time scale r of the phase fields: A certain time is necessary to establish local equilibrium against the 
"pull" of the advancing isotherms. 

A more surprising result is the rotation of the trijunction as a whole, which appears as soon as the symmetry between 
the two solid phases is broken. Here, we have shown only results at the eutectic composition for an asymmetric phase 
diagram; however, the same effect is recovered in the symmetric model alloy at off-eutectic composition, i.e., for 
unequal volume fractions of the two solids. The rotation does not seem to vanish when the interface width tends to 
zero. Although we cannot, at present, give a detailed explanation of this phenomenon, two speculative ideas about 
its origin might be pursued in the future. The first is that this effect could result from the local equilibrium of curved 
diffuse interfaces: Two different curvatures on the two sides of the trijunction could lead to some "torque" on the scale 
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FIG. 16: Amplitude ol saturated oscillatory limit cycles versus lamellar spacing, for parameters corresponding to the first 
(circles) and third (squares) lowest resolution in Fig. 115b (X/W — 64 and 128 there, respectively). Inset: same data, where the 
squares have been shifted along the A axis by 0.07. 



of the diffuse trijunction. The second is that the solute diffusion field in the liquid "wedge" close to the trijunction 
yields a torque term via its coupling to the diffuse interfaces, which would be a purely non-equilibrium effect. 

The origin of corrections to the classic sharp-interface model that do not vanish for sharp interfaces can be qual- 
itatively understood. Contrary to the situation around a dendrite or cell tip, where the lower cutoff scale for the 
diffusion field is set by the local radius of curvature, a trijunction point is a singular point in the sharp-interface 
solution. Indeed, when the flux lines of the chemical components are drawn as sketched in Fig. 1171 it can be seen that 
the curvature of these lines diverges as the trijunction point is approached. Therefore, the sharp-interface solution is, 
in some sense, inconsistent, because it is implicitly based on an assumption of scale separation between interface and 
"macroscopic" scales. Physical interfaces are smooth, and hence the local diffusion field around a trijunction varies on 
the same scale as the interface thickness. Therefore, the appearance of effects like the trijunction rotation is, in itself, 
not surprising. However, understanding them more quantitatively remains a challenging problem. More simulations 
are needed to study in detail the dependence of this phenomenon on the various parameters such as growth speed, 
contact angles, phase diagram and phase fractions. 

Recently, the dynamics of trijunctions was found to have a noticeable influence on the stability of lamellar arrays 
at small spacings |28|. More precisely, the stability threshold for lamella elimination was found to be distinctly 
smaller than the minimum undercooling spacing, contrary to previous predictions 0, |3^ . The reason is that the 
assumption used in these predictions that trijunctions move in the direction normal to the large-scale envelope of the 
composite front (Cahn's rule) was not satisfied. It should be stressed that this effect is not necessarily related to the 
trijunction behavior found here, since Cahn's rule does not make any assumptions on the local configuration around 
a trijunction. 



V. CONCLUSIONS 



We have presented and tested in detail a new phase-field model for two-phase solidification. The specific design of 
our free-energy functional has allowed us to clearly separate, at a local level, the dynamics of interfaces from those of 
trijunctions, and to eliminate all thin-interface corrections to the desired free-boundary dynamics for each solid-liquid 
interface. Both properties have been checked by numerical simulations: By the first we mean that interfaces between 
any two phases are completely free of the third phase far enough from the trijunctions, and this is indeed verified both 
at equilibrium and for moving interfaces (see Fig. [71) . As for the second, it implies that simulations at fixed physical 
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FIG. 17: Sketch of the diffusion field around a trijunction in the sharp-interface picture. Lines with arrows indicate the flux of 
solute around the trijunction. Sufficiently close to the trijunction, the concentration field obeys the Laplace equation, which is 
known to have a singularity in a wedge geometry. 



parameters become independent of the interface thickness, provided that it is not too large. To our knowledge, it is 
the first time that the quantitative convergence to a sharp-interface solution as the interface thickness is decreased, as 
shown in Fig. 1101 has been explicitly demonstrated for a phase-field model of solidification with more than one solid 
phase. 

For generic, non-symmetric situations (unequal volume fractions or physical properties of the two solid-liquid equilib- 
ria), small differences exist between the converged phase-field results and boundary-integral simulations. Furthermore, 
while the convergence is very satisfying for steady-state solutions, time-dependent morphologies such as oscillatory 
limit cycles still exhibit thin-interface effects. Whereas these prevent us, for the moment, from achieving calculations 
that are completely independent of the computational parameters, we have demonstrated that the resulting errors 
are very small, except possibly in the vicinity of bifurcation points. Since our model is not only very accurate, but 
also efficient, it is a suitable tool for the investigation of pattern formation in two-phase solidification. In particular, 
it can be applied to simulate three-dimensional structures and to study their morphological stability |43| . 

Because the individual interface dynamics are perfectly distinct and controlled in our model, the only possible source 
of the observed differences to sharp-interface models lies in the trijunction region. Indeed, our detailed numerical 
investigation has revealed that the behavior of the diffuse trijunctions in the phase-field model differs slightly from the 
assumptions usually made in sharp-interface theories, as already discussed in Sec. IIVFI Therefore, a quite surprising 
outcome of this study is that the calculation of two-phase microstructures is far more sensitive to the structure and 
dynamics of the trijunction points than one might have expected. Note that these findings have been possible only 
thanks to the specific design of our model, which has allowed us to eliminate all thin-interface corrections from the 
dynamics of the solid-liquid interfaces, and to separate them from those of the trijunctions. Therefore, our model 
is a privileged tool to explore in detail the physics of diffuse trijunctions, an interesting but challenging task. One 
conclusion of our study is that, apart from its undeniable fundamental interest, this more detailed knowledge is also 
a necessary ingredient for a future quantitative modeling of solidification with an arbitrary number of phases. 

Let us conclude by sketching how our model can be extended to treat more general situations. Here, we have limited 
ourselves to eutectic diagrams with straight liquidus and solidus lines and constant concentration jumps; however, the 
structure of our model makes it possible to simulate any phase diagram. In particular, phase diagrams with curved 
coexistence lines as well as peritectic phase diagrams can be simulated without difficulty by choosing appropriate 
coefficients A,{T) and Bi{T). 

Next, because of our specific choice for the coupling between composition and phase fields, the model lacks the 
effect of the local interface curvature on the solute rejection. This (negligibly small) effect could be recovered by the 
introduction of more complicated interpolation functions, as shown for single-phase solidification in Ref. 

Furthermore, the asymptotic analysis and the antitrapping current in its present form are only applicable if the 
two solid-liquid surface tensions are equal. In order to keep the benefits of the model for arbitrary surface tensions, 
the asymptotic analysis should be repeated (in part numerically) along the lines of Ref. as discussed in Sec. IIII HI 
this is a straightforward albeit tedious calculation. 

We expect that crystalline anisotropy can be incorporated in the model in the standard way by making the interface 
thickness depend on the local orientation of the interfaces. However, care has to be taken in oder to ensure that the 
anisotropic model still has interfaces that are free of third-phase adsorption. Finally, some considerations on how to 
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FIG. 18: Sketch of the phase-field state space. The representative point of the system is confined to the plane pi +p2 +p3 — 1, 
whose intersection with the axes defines the vertices of the shadowed triangle, the so-called Gibbs simplex. 

generalize the model to more than three phases are presented in appendix lUl 
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APPENDIX A: TIME CONSTANTS 

Here, we show that a linear relation between the time derivatives of the phase fields and the driving forces does not 
permit to choose different relaxation times for different interfaces and still guarantee that each interface is free of any 
third phase. 

A given configuration of the phase fields can be visualized by a representative point in the three-dimensional vector 
space spanned by the axes pi, P2 and (see figure [T5)l . The constraint P1+P2+P3 = 1 restricts the allowed positions 
to a two-dimensional subspace, the plane that contains the points (1, 0, 0), (0, 1, 0), and (0, 0, 1). These points are also 
the vertices of an equilateral triangle, the Gibbs simplex (shadowed area). Only points inside the triangle (shadowed 
area) are meaningful if each phase field is to be interpreted as a positive volume fraction. 

The most general form for a linear relationship between the time derivatives of the phase fields and the driving 
forces SF/Spi is 

t-Sr«£. (AD 

j=l •^■1 

where Fy can be seen as a linear map between the vector of the driving forces, 6F/dpi, and the vector of the time 
derivatives, dpi/dt. Since the representative point must remain in the plane which contains the Gibbs simplex, the 
vector of the time derivatives must also lie in it, and therefore be normal to the vector (1,1,1). In contrast, there are 
no explicit constraints on the vector of the driving forces, which can have any direction. Therefore, the matrix F^ 
must have the eigenvector (1, 1, 1) with eigenvalue zero, or, in other words, be a projection on the allowed plane. 

In addition, the free-energy density F is designed so that on purely binary interfaces there is no driving force for 
the formation of the third phase. For example, on the 1-2-interface, the driving force along the p3-direction is zero; 
therefore, the vector 5F/Spi is parallel to (1,-1,0). Since we do not want the resulting time evolution to introduce 
the third phase, the vector of the time derivatives must also be parallel to (1,-1,0). This implies that this vector 
must be an eigenvector of the matrix Fy . The corresponding eigenvalue is the relaxation rate for the 1-2-interface. 
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The same reasoning can be repeated for the 1-3-interface, with the associated eigenvector (1,0,-1). Since the 
two vectors are not orthogonal, the only possibility that makes both of them eigenvectors is that the whole subspace 
corresponding to the plane of the Gibbs simplex is an eigenspace with a single eigenvalue. Therefore, all the time 
constants associated with purely binary interfaces are the same, and the only valid choice for the matrix F is a multiple 
of the projector on the plane of the simplex, 

r,y =r(% (A2) 

where u — (1, 1, l)/\/3 is the unit vector normal to the simplex plane. Writing F = 1/{tH), the resulting equations 
of motion are 

This is equivalent to Eqs. (|3.4|) in conjuction with the prescription of Eq. H3.5|) . 

The only way to achieve different relaxation times for the different interfaces while keeping purely binary interfaces 
is therefore to introduce some nonlinearity in the relation between time derivatives and driving forces. We achieve 
this by letting the time constant r depend on the phase fields. Note that locally "anisotropic" dynamics in the 
configuration space could also be induced by letting all the coefficients depend on the phase fields, so that the 
directions of the eigenvectors and the associated eigenvalues depend on the location in the configuration space. 



APPENDIX B: SYMMETRY OF THE GRADIENT ENERGY COEFFICIENTS 

The gradient energy expression, Eq. (|3.11|) . satisfies the conditions for no third-phase adsorption, Eqs. (|3.9f) . We 
now consider its simplest generalization, 

/grad ^Y^^P"^ • ^P^^ (Bl) 

mn 

a quadratic form in the gradients of the phase fields, where the ^„in are dimensionless constants that we assume to be 
symmetric (and positive), ^mn = £,nm{> 0) Vm, n. We find that Eq. H3.9al) is satisfied if and only if the coefficients that 
multiply the Laplacians of the phase fields in each purely binary interface are all the same, so that this generalization 
does not provide any extra freedom for our purposes. 

To see this, compute the functional, constrained derivatives of a free-energy functional -Fgrad defined as the volume 
integral of the gradient terms alone, Fg^ad — Jy -f^'/grad- Making use of Eqs. H3.5|l and H3.6|l . we find 



rad 



Uk^^Pm ^^nn^^Prn ) • (B2) 



Pq+P/3+Pl = 1 V m ^ mn 



Spk 

The condition of flatness with respect to variations of pk on a i-j interface, Eq. I|3.9a|l . then yields 

- 2^iA; = Cjj" ^ 2^ifcj (B3) 

where we have taken into account that pk = and pj ^ I — pi . This can be written in a more symmetric way by 
adding S^kk to both sides: 

^kk — 2^ife + in ~ £.kk — "^ijk + ijj- (B4) 

Two equivalent conditions can be obtained for the other two interfaces, so that all three symmetric combinations 
imm — '^S.mn + ^nn must havc the Same value. 

Next, replace k with i in Eq. ljB2|) and, afterwards, apply it to a i-j interface {pk = 0). One finds 



Spi 



= Y (2^- + - SCy - i^k + ijk)y^P^ ■ (B5) 



Using Eq. 1B3|I to eliminate ^jk — S,ik, we find that the prefactor of the Laplacian can be rewritten as {K/2)[iii — 
2^y -I- ijj), i.e., one half of the combination that we saw above to be the same for all interfaces. Therefore, we can 
choose = Sij without loss of generality, and we are back to the minimal model with Eq. (|3.11l) . 

The multi-phase-field models based on the approach of Steinbach et al. [2^ use a different expression for the 
gradient terms. It combines gradients of the phase fields with the fields itself, and is hence not a simple quadratic 
form in the gradients. We have checked that this expression does not satisfy the condition of flatness; therefore, it is 
not a suitable alternative for our purposes. 
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APPENDIX C: MORE THAN THREE PHASES 



In principle, the multi-phase formalism where each phase field represents a local volume fraction can be extended 
to an arbitrary number of phases. However, the methods we have applied to ensure that interfaces are free of any 
adsorbed third phase need to be generalized with care. Let us consider a multi-well potential /mw constructed in 
the same manner as our triple- well potential /tW: i-e., by a superposition of individual double- well potentials for 
each phase field, /mw = J2i fnwiPi)- For fiywip) — bl"|l ~ (with v a positive real) the free energy of a point 
where n phases coexist (n-junction) is n^~'^^(n — 1)". For our model, v = 2, this expression has a single maximum 
at n = 3. Therefore, a quadrij unction already has a lower free energy than a trijunction, which means that the 
free-energy landscape /mw has undesirable local minima for more than three phases. These lower-energy n-junctions 
can always be "raised" by adding one obstacle term /obs — ^i'=iPi per rt-junction to treat. However, the number of 
such terms becomes rapidly prohibitive when n increases; in addition, such "steep" terms require a higher numerical 
resolution. Alternatively, for < < 1 the n-junction energy becomes monotonously increasing in n for all positive 
n. Therefore, the simplest choice ensuring the right topography for an arbitrary number of phases would he v = 1, a 
model considered in Ref. |3ll |: however, to our knowledge no thin-intcrface analysis has been performed to date for 
this inverse parabolic potential. 
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